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ABSTRACT 

The three-dimensional linearized vorticity transport 
equations for plane and pipe Poiseuille flow were studied 
using a highly generalized complex exponential form of 
solution in both space and time. The stability of these 
flows was examined using frames of reference which move with 
the fluid particles. 
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TABLE OF SYMBOLS 



All guantities below except those subscripted with d are in 
dimensionless form. The nondimens ionalization is described 
in Appendix A. 



D,D 2 , 


The partial derivatives with respect to y in 
cartesian coordinates or with respect to r in 
cylindrical coordinates. 


e 


Base of natural logarithms. Exponentiation is 
also denoted by exp ( ) 



e # e ,e Unit vectors along the x, r, and 0 axis in 
x r 0 





cylindrical coordinates. 


f ,g,h 


Components of the velocity vector potential 
defined in eguations (2-9a) and (E-1) . 


i 


+ V- 1, the imaginary unit. Also used as index 
in finite differencing mesh in Section III. 


i* j/k 


Unit vectors along the x, y, and z axis in 
cartesian coordinates. 


L 


Reference length. Radius of pipe in 

cylindrical coordinates and semiheight of 

channel in cartesian coordinates defined in 
Appendix A. 


n 


The number of interior points in the finite 

differencing mesh of Section III. Also used 

as imaginary part of $ in pipe flow. That is, 

/3 = /5 + in. 

R 
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p 


Pressure. 


Be 


Reynolds number based on mean velocity and 
channel semiheight or pipe radius. 


t 


Time. 


T,t 


Shorthand notation for commonly occuring 
groups of symbols defined in eguations (D-54) , 
(E-34) and (E-48) . 


U,V,V 


Components of the complex perturbation 

velocity defined in equation (2-9b) . 


u ,v ,w 
n n n 


Components of arbitrary vector function V . 

n 


U* / V* / H* 

n n n 


Components of V*. 

n 


V 


Velocity vector of the perturbation flow 

defined in Appendix A. 


V 


Velocity vector of unperturbed laminar flow 
defined in equations (2-3) and (2-5) . 


V 

avg 


Mean velocity defined in Appendix A. 


V 

n 


Arbitrary vector function in cylindrical 

coordinates with angular periodicity n. Used 
in Appendix G. 


V* 

n 


Vector function V rotated $> about the origin. 

n 


W 


Complex vector potential of perturbation 
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I u 



x,r,0 


velocity defined by equation (2-7) . 

Cylindrical coordinates defined in Figure 1-2. 


x / y / z 


Cartesian coordinates defined in Figure 1-1. 


V ,p ,t 
d d d 


Dimensional variables of velocity, pressure 
and time used in Equation (1-1) . 


a 


a + ia Complex wave number of the 

R I 

perturbation in the x direction. 


P 


B + i $ Complex wave number of the 

R I 

perturbation in the z or 0 direction. 


r 


y + i Y Complex frequency of the 

R I 

perturbation . 


f 

r ,r , r 

x y z 


The vorticity transport equation (D-55) 

expressed in abbreviated notation as defined 
by equations (B-1) and (B— 11). 

The components of T in cartesian coordinates 
defined in equation (B-1). 


r,r , r 

x r Q 


The components of T in cylindrical coordinates 
defined in equation (B-11). 


V 


Kinematic viscosity. 




Components of the complex perturbation 

vorticity defined in equation (2-9c) . 
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Density in equation ( A— 1 > - 

Angle of rotation of arbitrary point about the 
origin in cylindrical coordinates as in figure 

G-1 . 

Vorticity vector of the perturbation flow 
defined in eguation (A-8) . 

Vorticity vector of unperturbed laminar flow 

defined in equations (2-4) and (2-6) . 

Linear vector operator (nabla) . Used for 
gradient (v) , divergence (v*) and curl (vx) . 

Sign of vector cross multiplication. 

Brackets enclosing a matrix. 

Brackets enclosing a column vector. 



I. background 

The problem of solving the Navier-St okes equation for 
the onset of instabilities in laminar flow for simple 
geometries has been studied for many years. The problem is 
reduced to its simplest form by linearizing the vorticity 
transport equation, assuming sinusoidal perturbations in 
space and complex exponential perturbations in time, and 
reducing the eguations to their two-dimensional form. The 
solutions to th is simplified problem are too stable to agree 
with experimentally obtained values of the critical Reynolds 
number. 

To compare Reynolds numbers from dissimilar flows, the 
hydraulic diameter and mean flow are generally used as the 
characteristic length and velocity. Using these criteria, 
the critical Reynolds number for a pipe is 2600 [Schlichting 
1968] and for a rectangular channel with a width to height 
ratio of 8:1 is also 2600 [Kao and Park 1969], The analysis 
in this paper uses the channel semiheight and the pipe 
radius for analytic convenience (with the mean velocity) . 
Based on these characteristic lengths the experimentally 
obtained critical Reynolds numbers become 1300 for a pips 
and 730 for the 8:1 channel. According to Kao and Park the 
critical Reynolds number for plane Poiseuille flow will be 
greater than that for a channel of finite width. 

The possible destabilizing effect of three-dimensional 
flow in the case of plane Poiseuille flow was ruled out by 
Sguire's theorem [Sguire 1933] which proved that 
three-dimensional perturbations become more unstable as they 
are rotated into the plane of the two-dimensional flow. 
Generalizing the problem by allowing exponential growth of 
the perturbations in space as well as time greatly adds to 
the complexity, and the solution of the general problem has 
not been previously approached. 
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Much work has been done for the case of pipe flow. 
Flows with various angular wave numbers, and with sinusoidal 
spatial perturbations were shown to be stable [Salwen and 
Grosch 1572], perturbations with exponential growth in 
space, but with strictly sinusoidal time variation were 
studied and found stable [Garg and Rouleau 1971] and the 
combination of exponential growth in space and in time was 
studied using power series analysis [Gill 1973] and all were 
found to be stable. 

Because of the general agreement that pipe flow is 
stable to infinitesimal disturbances, two other approaches 
to the problem of accounting for the critical Reynolds 
number have been used. First, it has been assumed that 
perturbations which are stable while they are infinitesimal, 
become unstable when they are finite. Studies have 
concluded that finite disturbances are unstable for pipe 
flow [Davy and Drazin 1968]. Similar theoretical 
examination of plane flow has reached the same conclusion, 
that is, finite perturbations are less stable than 
infinitesimal perturbations [ Hclntire and Lin 1971]. 
Second, it has been postulated that the origin of the 
critical instabilities occurs in the entrance region of the 
pipe (or channel) and both theoretical and experimental 
studies have been done to demonstrate this [Huang and Chen 
1973] [Leite 1956]. 

These two explanations for the observed instability of 
Poiseuille flows have been postulated as a result of the 
inability of the linear theory to account for the 
experimental facts. Unfortunately, a totally general 
solution to the linear problem has never been accomplished. 
In this paper such a solution is presented. Perturbations 
which have a fully complex exponential form and which are 
fully three-dimensional are introduced and the resulting 
equations are solved using finite-difference methods. 
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Results for channel flow indicate that by varying the real 
part of the spatial wave number, the corresponding critical 
Reynolds number can be lowered. This is promising because 
it might help resolve discrepancies between theory and 
experiment. Further study is necessary to clarify these 
results. 
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II. theory 

A. EASIC EQOATIONS 

The unperturbed laminar flow of an incompressible fluid 
with constant viscosity is governed by the vorticity 
transport and continuity eguations. These involve the 

velocity and vorticity vector functions V and Cl which are 

well known for the two simple geometries being considered 
(see eguations (2-3) through (2-6) ) . If the flow is 
slightly perturbed and the perturbation velocity and 

vorticity are denoted by v and co respectively, the 

perturbation continuity equation and the linearized 
vorticity transport equation are easily derived (see 
Appendix A) . These equations, in their nondimensional form, 
may be written 

v*v = 0 (2-1) 



and 



1 v 2 u> + o>*vV - V«vw + v®v£l 
Ke 

- &®vv - 9ai/3t = 0 (2-2) 

where w = vxv. For the two cases under consideration, plane 
Poiseuille flow and pipe flow, the coordinate axis 
orientation is given in figures 2-1 and 2-2 with the 

velocity and vorticity distributions V and Cl. 
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For Plane Poiseuille Flow 



V = i3(1-y2) = iU(y) 


(2-3) 


ft( y) = vxV (y) = k3y 


(2-4) 




Figure 2-1 Plane Poiseuille Flow 



For Pipe Flew 

V (r) = e U (r) = e 2(1-r2) (2-5) 

x x 

i2(r) = vxV (r) = e 4r (2-6) 

e 




Figure 2-2 Cylindrical Poiseuille Flow 
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The vorticity transport equation in three dimensions is 
actually three separate equations. These three equations 
are not independent (as shown in Appendix B) and represent 
only two fully independent conditions. The continuity 
equation is a scalar equation and the relation between 
vorticity and velocity is a vector equation. Thus there are 
six independent equations in six unknowns. The unknowns are 
the three components of the perturbation velocity and the 
three components of the perturbation vorticity. 

These can be reduced to three equations in three 

unknowns by replacing cu in equations (2-1) and (2-2) by vxv„ 
A further simplification can be made by introducing the 

velocity vector potential, W, where 

v = vxW . (2-7) 

The continuity equation is satisfied identically because 

v«v = v (vxW) = 0 (2-8) 

by a well known vector identity. If v is now replaced by 

vxW, the result is two equations (the two independent 
conditions cf the linearized vector vorticity transport 
eguation) in three unknowns. The unknowns are the three 
components cf the velocity vector potential. From this 

result it may be deduced that the components of W are 

redundant. This is indeed the case and, as shown in 
Appendix C, one cf them may be set arbitrarily to zero. 

A useful way to take advantage of the linearity of 
equation (2-2) is by seeking solutions which are complex. 
Both the real and imaginary parts of any solution obtained 
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will, by themselves, be solutions. By choosing solutions of 
an exponential form, as in equations (2-9), the system of 
partial differential equations is reduced to ordinary 
differential equations. The desired form for solutions in 
cartesian coordinates is 



W (x , y, z,t) = £ if (y) + 


X 

jg (y) + kh (y) ]e 


(2-9a) 


v (x,y, z,t) = £ iu (y) + 


X 

jv (y) + kw (y) ]e 


(2-9b) 


"(x,y,z,t) = £ i £(y) + 


X 

j77(y) +£k (y) ]e 


(2-9c) 


where 






X = ax+/3z+yt 




(2-10) 


For the cylindrical case 


, the same equations. 


(2-9) and 


(2-10), apply with x,y,z 


replaced by x,r,0. 


The wave 


numbers a, /3, and y are, in 


general, complex. The 


functions 



f, g, and h are defined by equation (2-9a) . They are the 
part of the components of W which contain the y (or r) 
dependence and, in general, are complex. Similarly, u, v and 
w are part of the components of u. And £, rj and £ are part 
of each component of w, all complex functions of y (or r) . 

B. CARTESIAN COORDINATES 

The plane flow case for cartesian coordinates is 
considered first. Using equations (2-9) to substitute into 
equation (A-13) , as described in Appendix D, results in an 

equation in terms of W. 
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-1 vxvxvxvxW + vxVxvxvxW - vxiixvxW 

He 



- vxvx9 W = 0 



(D-1) 



Appendix E then develops the form of the coefficients of 
equation (D-1) as matrices. The final equation is 

Cd^ (D3f') fD 2 g 

e:§ + £m 3 ] rM + ([H 2 ^ CN a ]) fei 



+ cr« l 3+rt ]) jg|r + ([ m q ]+x£ n q ]) 



= o 



(D-55) 



Where D, D 2 , ... are the differential operators with respect 
to y. The matrix coefficients are given in equations (D-46) 
through (D-54) . This vector equation may be expressed in 
abbreviated form by 



r = 



= o 



(B-1) 



as discussed in Appendix B. 
three separate equations 

r = 0 

x 



Equation (3-1) is actually 



(B-2a) 



r = 0 



(B-2b) 



r = o 



(B-2c) 



Examination of 
(D-55) shows that 



the matrix coefficients of equation 
the highest order derivatives of the 
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\ 




components f, g, and h, of W in each of the three equations 
is as fellows 



eguation 



r 

x 



r 

y 



r 

z 



highest 
order of f 

4 

3 

2 



highest 
oraer of g 

3 

2 

3 



highes t 
order of h 

2 

3 

4 



TABLE 2-1 



As demonstrated in Appendix C, one of the components of 

W may be set identically to zero, but, the choice of h(y)=0 
leads to a degeneracy in the eguations for the classical 
case of /3 = 0. Since it is desirable to have the 

three-dimensional case reducible to the plane flow, 
two-dimensicnal case, it is stipulated that either f or g be 
set to zero. Appendix F develops the form of the boundary 
conditions for the three cases of setting f, g, or h to 
zero. For g{y)=0, the boundary condition 

f(±1) = h (± 1) (F-6a) 

is more complicated than for the case of f(y)=0. It is 
therefore desirable to set f(y) to zero. For f ( y ) = 0 , the 
boundary conditions are 

g (± 1) = 0 (F-4a) 

h(±1) = 0 (F-4b) 
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Dh (± 1) = 0 . 



(F-4c) 



From table 2-1, derivatives of f (y) occur to the fourth 
order. Eguations (F-4b) and (F-4c) give four boundary 
conditions for f (y) . Table 2-1 also shows that g(y) occurs 



to the third order in eguations F and T , while there are 

x z 



only two boundary conditions in eguation (F-4a). It is 



necessary tc reduce the order of g in F and T . This can 

x z 



be done by taking a linear combination of T and T . 

x z 



Examination of the matrix of eguation (D-47) shows that the 



coefficient of D 3 g in eguation T is a. In eguation T the 

x Re z 

coefficient of D 3 g is P_. Thus, the following combination 

Re 

of eguations F and F will eliminate the terms in D 3 g. 
x z 



-PF + otT = 0 (2-11) 

x z 



Appendix B shows that this linear combination of eguations 



r and T , and the condition 
x z 



F = 0 (2-12) 

y 

are sufficient to satisfy the vorticity transport eguation 
written as in eguation (B-1), under the condition that 

a 2 + p z * 0 . (B-9) 



The resulting eguations to be 
(2-12) in the unknowns g (y) and h (y) . 



solved are (2-11) and 
Performing the actual 
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calculations indicated by equation (2-11) reveals an 
interesting and extremely convenient fact. Not only the 
third order derivatives of g disappear from equation (2-11) , 
but g and all its derivatives cancel out of the equation. 
Setting f (y) to zero, equation (2-11) only contains h (y) . 
Thus, the eguaticns are uncoupled and (2-11) may be solved 
for h (y) alone. With h (y) known, equation (2-12) may then 
be solved for g (y ) . Writing out equation (2-11) yields 

-_a_D*h + (aT-^a_ (a 2 + yQ 2 ) )D 2 h + (3a 2 + a(a 2 +/3 2 ) T) h 

+ r (aD 2 h + a(a 2 +/3 2 )h) = 0 (2-13) 

where 



T = aU - 1 (a 2 + /3 2 ) = 3 a ( 1-y 2 ) - 1 (a 2 +fi 2 ) (D-54) 

Be H 2 Be 

The boundary conditions for equation (2-13) are 

h (± 1 ) = 0 (F-4b) 

Dh (± 1) = 0 . (F-4c) 

This homogeneous differential equation is solvable as an 
eigenvalue problem. This technique is discussed in the 
section on numerical methods. The associated problem is to 
solve equation (2-12) for g(y). The coefficients of (2-12) 
may be obtained directly from equations (D-48) through 
(D-54) . This inhomogeneous equation may be solved by 
various integration techniques to obtain the associated 
function of g (y) for each of the eigenvalues obtained from 
the solution of equation (2-13). Solving equation (2-12) 
would be cf interest if the functional shape of the 
perturbation velocities and vorticities were desired, but 
the problem was net dealt with in this paper. 
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c. 



CYLINDRICAL COORDINATES 



The equations for pipe flow, in cylindrical coordinates, 
are developed in Appendix E. Equation (D-1) is still valid 
and in abbreviated form may be expressed by 



s \ 



r =< 



r > 

r 






= o 



(B-11) 



as discussed in Appendix B. This vector equation is 
actually three separate equations. 

T = 0 (B- 1 2a) 

x 



r = 0 (B- 1 2b) 

r 



e 



o 



(B- 1 2c) 



The final matrix form of equation (B-11) is given in 
equation (E-55) with the coefficient of equations (E-40) 
through (E-48) . Examination of these coefficients shows 
that the highest order of derivatives of the components f, 

g, and h, of H in each of the three equations is as follows 
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highest highest highest 

equation order of f order of g order of h 



‘0 



TABLE 2-2 



As in the cartesian case, the order of g(r) in equations T 



and T is too high for the boundary conditions. To 
6 

the order of derivatives of g (r) in T and T a 

x 0 

combination is taken. Checking the coefficients of 



both eguations it is found that the combination 



-§r + ar = o 

r x 0 



x 

reduce 
linear 
D 3 g in 

(2-14) 



eliminates the third order derivative of g (r) . This 
combination of eguations does not, in general, uncouple the 
problem by eliminating one of the components of the velocity 
vector potential from the equation. The resulting matrix 

equation, with all three components of W (one of them may 



later be set to zero) can be represented in the same form as 
eguation (D-55) . 





([M i ] + /[N i ]) 



Df') 



f D2 g) 

(CM a ]*y£H 2 ])jD«g| 

+ ([ M o ]+rl N 0 3) jl] 



= 0 



(D-55) 
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but now the matrix coefficients are 2X3 matrices, where the 

top row is the combination of eguations given in (2-14) and 

the bottom row is eguation T . The coefficients are as 

r 

follows . 



[H ] 

4 



' 8 

rHe 

0 



0 

0 



- a 
TTe 



(2-15) 



[« 3 ] 



28 


0 


- 2a 


File 




FRe 


a 

He 


0 


rfe 



(2-16) 



[M 2 ] = 



/St- 8 -/St -4 a/3 
r he r 3 Re r F^TSe 

a -t 

rHe He 



aT + 3a -at 

r^lTe He 

2/S 

r^He 



(2-17) 






- /St- 383+ 8 + a 2/S 

r 2 rTKe Fie r^Re 



-aT- a 
r^He 



4a8 aT+ 3a8 2 - 3a - a 3 

r 3 Ee r r 3 Ke r r We FRe 



8 2 - as 
r 3 Re rHe 



-^T- _8_ 
r r^Ke 



£M o ] 



-£tT + 483 2a8T-4a8 -2a8t 

r r b Re r 2 F’HTe r 2 Re 
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( 2 - 22 ) 



£N o ] = 



-fit 

r 



2 afi at- J 
r -z r ? 



0 



t 



-fi 

r? 



where 



T = aU 




2 a(1-r 2 ) - 1 (a 2 +fif) 
He r* 



and 




(E-48) 



(E-34) 



Since the parameter fi, which determines the angular 
periodicity, only has discrete values as determined by the 
boundary conditions in Appendix G, 

fi = ni, n=0,1,2,... (G-4) 

each value of fi is considered independently. For fi=0, it is 
found from the coefficients (2-15) through (2-22) that the 
eguations uncouple, even more completely than in the 
cartesian coordinate case. In the first equation 
represented by these coefficients (the top row of the 
matrices), both g (r) and f(r) drop out. The coefficients 
remaining represent a fourth-order ordinary homogeneous 
differential equation in h (r) . The boundary conditions for 
this problem are derived in Appendices F and G. Taking f (r) 
to be identically zero, the boundary conditions at the wall 
are 



h ( 1 ) = 0 (F- 15b) 

Dh (1 ) = 0 (F- 15c) 

The boundary condition at r=0 from (G-40) is that all even 
derivatives of h are zero. In the second equation h(r) 
drops out because all coefficients contain fi. With f set 
identically to zero, this equation becomes a homogeneous. 
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second order differential equation in g (r) . The boundary 
condition at r= 1 is 

g ( 1) = 0 (F- 1 5a) 

The boundary condition at r=0 from (G-40) is that g(0) and 
all even derivatives of g are zero at the origin. 

The consequence of the complete separation of the 
components g and h into two separate homogeneous 

differential equations is that there are two sets of 
eigenvalues which are derived independently of one another. 
The velocity in terms of f, g, and h is derived in Appendix 
F. With >3=0 and f(r)=0, equations (F-13) become 

u (r) = JD (rh (r) ) (2-23a) 

r 

v (r) = ah (r) (2-23b) 

w (r) = ag (r) (2-23c) 

Thus the axial and radial components are related to the 
function h (r) , while the angular velocity is completely 
independent . 

Both of these homogeneous differential equations can be 
solved as eigenvalue problems. This technique is discussed 
in the section cn numerical methods. From the eigenvectors 
of f and g the velocity and vorticity perturbation forms may 
be found. This was not done in this paper. 

If n=1, that is, .if /3 = /? + i, then the two equations 

R 

represented by eguations (2-15) to (2-22) do not uncouple 

and must be solved simultaneously, with two of the 
components of the velocity vector potential as the unknowns. 
From Appendix F it may be observed that the most complicated 
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boundary condition at the wall 

£f (1) = ah (1) (F- 17a) 

arises from setting g (r) egual to zero. From Appendix G it 
is also found that one of the boundary conditions for n=1 at 
r=0 is 

g (0) = ih { 0) (G-41) 

These inhomogeneous boundary conditions may be avoided by 
the choice of setting h(r)=0. 

The result is that the first and second columns of the 
matrices of eguations (2-15) through (2-22) give the 
coefficients of f and g for each of the two eguations (the 
top and bottom rows of the matrices) . The boundary 
conditions at the walls are 

g ( 1) = 0 (F-1 9a) 

f(1) = 0 (F- 19b) 

Df { 1) = 0 (F- 19c) 

The boundary conditions at r=0 are that f(0) and all even 
derivatives of f are zero, and that g(0) and all odd 
derivatives of g are zero. 

Thus, the problem for n= 1 becomes two simultaneous 
homogeneous differential eguations which must be solved 
together. This problem is discussed in the section on 
numerical methods. 

B. STAEILI1I C RITEBION 

Previous investigations into the guestion of the 
stability of Poiseuille flows have used various criteria for 
determining the stability of the flow from the solutions 
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obtained. In particular, the real part of the eigenvalue y 

R 

is usually used when there is no real-exponential spatial 
variation [Salwen and Grosch 1972]. Where exponential 
growth in space has been a part of the problem [Garg and 
Rculeau 1971] the real part of the spatial wave number has 
been taken to give the instability. In other cases, the 
phase velocity has teen used to determine the stability. 
Clearly, a measure of the stability of the flow must take 
into account both the exponential rate of growth with time 
and with space. 

Since it is the stability of the perturbed fluid 
particles that is of concern, it is natural to consider how 
the rate of growth or decay of the perturbations is seen 
from a coordinate frame moving with a particular fluid 
particle. For cartesian coordinates, the fluid particles 
can have velocities ranging from 0 to 1.5 depending on their 
distance from the walls, that is, depending on the value of 
y. The distribution of fluid velocities is the Poiseuille 
distribution of eguation (2-3) . 

U(y) = 3 ( 1-y 2 ) (2-24) 

7 

These moving axes will have a velocity with the mean flow in 
the x direction. Let x', y, z, t be the coordinates and a, 
P, y* the complex wave numbers with respect to the moving 
axes. The form of the perturbation vector potential for a 
given eigenvalue obtained as a solution is 

W = (39 (y) + kh (y) ) exp (ax + Pt. + yt) . (2-25) 

The complex frequency y' seen from this moving reference 
frame will be different than y from the fixed frame. To 
establish the relation between y' and y, the perturbation is 
written in the moving frame and then transformed into the 
form which corresponds to the fixed frame. 
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H 



(jg + kh)exp(ax' + $z + y't) 



(jg + kh) exp[ a (x-Ut ) + /3 z + y* t ] 



(jg + kh)exp[ax + @z + (x'-aU)t] 



(jg + kh)exp[ax + /3z + /t] 



(2-26) 



Thus 



/' - a U = y 



(2-27) 



Solving for y % and splitting into real and imaginary parts 
gives 



If y % is positive, zero, or negative, the perturbation is 
R 

said to be unstable, neutral, or stable, respectively. 

Hence, the value of y' is taken to be a measure of the 

R 

stability of each eigenvalue obtained. 




(2-28) 




(2-29) 
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III. NUMERICAL METHODS 



A. CARTESIAN COORDINATES 



The equation to be solved is the linear combination of 
the first and third components of the vorticity transport 
equation which uncouples the dependence of h (y) and g (y) . 
This is developed in Section II and results in equation 
(2-13). 

-q_D*h + (aT- a_ (a 2 +£ 2 )) D 2 h + (3a 2 + a(a 2 +# 2 ) t) h 

+ y(aD 2 h + a (a z +p 2 ) h) = 0 (2-13) 

where T is defined by equation (D-54) . 

T = au - 1 (a 2 +£ 2 ) = 3<*(1-y2) - 1 (a 2 +£ 2 ) (D-54) 

He 2 He 

The boundary conditions are derived in Appendix F. 

h (± 1) = 0 (F-4b) 



Dh (± 1) = 0 . 



(F-4c) 



Various methods are applicable to the solution of this 
homogeneous fourth-order equation. The method of solution 
applied is the method of finite differences, approximating 
the function h (y) by n discrete evenly spaced unknowns, and 
solving for the eigenvalues and eigenvectors of the system 
of equations thus generated. 

One of the variables, a, P or y must be chosen to be the 
unknown eigenvalue. Since y is linear while a and $ occur 
to higher powers, the choice of y produces a more standard 
problem. Thus, y is the unknown and values must be assigned 
to a and P. Equation (2-13) has already been written in the 
appropriate form for this method of solution, with the 
factors of y separated from the rest of the equation. In 
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an abbreviated, but general, form, equation (2-13) can be 
written 

C D 4 h + C D 2 h + C h + y(C'D 2 h + C ' h) = 0 (3-1) 

* 2 0 2 0 

where the coefficients (C's) are, in general, functions of 
a , P, and y . 
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Figure 3-1 Finite Difference Mesh 

The range of y, from +1 to -1 across the channel height, 
is divided into a one-dimensional computational mesh of even 
spacings. As shown in Figure 3-1, there are n interior 
points, n+ 1 divisions and n+2 total mesh points, including 
the boundaries. The grid spacing is Sy. 

If the values for the function h (y) at the grid points 
i-2, i-1, i+1 and i+2 are expanded as Taylor series about 
th 

the i grid point, the resulting set of simultaneous 
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equations 


may be solved for the second- 


order central 


difference 


approximations of the derivatives 


. These are as 


follows . 






Dh = 
i 


(h. -h. )/2(8y) 

1+1 i-l 


(3-2a) 


D 2 h = 
i 


(h. -2h +h J/(Sy)2 

1+1 l 1-1 


(3- 2b) 


D 3 h = 
i 


(h. -2(h. -h. )-h. )/2(8y)3 

1+2 1+1 1-1 1-2 


(3-2c) 


D 4 h = 
i 


(h -4h +6h -4h +h )/(Sy) 4 

i + 2 i+1 i i-1 i-2 


(3- 2d) 



The error cf equations (3-2) is of the order of magnitude 
(Sy) 2 . 

The values of h at the boundaries h(±1) are known to be 

zero by the boundary condition of equation (F-4b) . Thus, 

there are n unkncwns--the values of h (y) at the n interior 

points — denoted h , h , h,...,h. By substituting 

12 i n 

equations (3-2) for the derivatives of h in equation (2-13) , 

the fourth order differential equation becomes a set of n 
linear, algebraic, difference equations which form a banded 
matrix with a band width of five. Each equation, written 
th 

for the i mesh point, includes the values of h at the mesh 
points i-2, i-1, i, i+1 , and i+2. There are four special 
cases which must be handled using the boundary conditions. 
The equations written for the points 2 and n- 1 include the 
boundaries. Since the values for h at the boundaries is 
known to be zero, those terms drop out of those equations. 
The equations written for points 1 and n contain terms for 
the boundaries and terms for virtual points outside the 
boundaries, which, in accordance with the labeling scheme. 
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Figure 3-2 Finite Difference Matrix 



may be called the -1 and n+2 values of h (see Figure 3-2). 
The second boundary condition, applied at the zeroth 
(boundary) point, yields the condition 




Thus, in the difference equation written at the first 
interior point, the term for the boundary is dropped and the 
coefficient of the virtual point is combined with the 
coefficiert at the first interior point. The equation 

th 

written for the n point is treated similarly. 



The matrix thus formed is broken into two matrices, one 
containing the coefficients of the eigenvalue y, and the 
other containing the remaining terms. The resulting matrix 
eguaticn may be written 

[X]jvj - rm{v] = 0 (3-4) 

Matrix [Y] will be a three-banded matrix because y only 
occurs as a factor of D 2 h and h. £vj is the eigenvector of 
values of the function h at the mesh points. Figure 3-3 
shows the structure of the matrices. 
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Figure 3-3 Matrix Form of Equation 
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The two subprograms, FUNCTION CHM1E1 and FUNCTION CHN2E1 
supply the coefficients of these matrices. The input 
parameters, K and Y are, respectively, the relative position 
in the band for each row, and the value of y for the 

difference equation. The subprograms contain function 

statements which produce the coefficients obtained from 
equation (2-13). The parameter K determines how those 
coefficients are combined to yield the coefficients of the 
points i-2, i- 1 , i, i+1, and i+2, which are determined using 
the central difference equations, (3-2) . Because the 
standard eigenvalue equation has a minus sign, as in 
equation (3-4), the signs of the coefficients of CHM1E1 are 
reversed from the sign of the terms in equation (2-13). 

The eigenvectors obtained from this problem have either 

even or odd symmetry, that is, they have symmetry about y=0 

such that either 

h (y) = h (-y) (3-7a) 

or 

h (y) = -h(-y) . (3-7b) 

Using this fact, the problem size can be cut in half by 

solving for the values of h across only half of the 

i 

channel, once with boundary conditions at y=0 corresponding 

to even symmetry and a second time with boundary conditions 
corresponding to odd symmetry. This gives the full set of 
eigenvalues and vectors. The boundary conditions for the 
mesh points near the center of the channel are easily 
derived from the conditions in equation (3-7) . 

Subroutine MSET sets up the coefficients of either 
matrix X or matrix Y. The variable MODE is used to control 
how the matrices are set up. If MODE=0, the matrices for 
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the full channel result. M0DE=1 or M0DE=2 will result in 
matrices set up for the half channel with boundary 
conditions for either odd or even symmetry. For the 
half-channel solution, the total number of divisions in the 
channel is always even so that the center of the channel, 
y=0, falls between two mesh points. The center of the 
channel could fall directly on a mesh point but this would 
result in the inconvenience of having one more eigenvalue in 
the solution for even vectors than in the solution for odd 
vectors. 

At the time of writing there were no computer routines 
available tc solve an eigenvalue problem with complex 
matrices in the form of equation (3-4) . Therefore, the 
problem is converted to the more conventional form 

([Z] - y[i J) V =0 

by inverting matrix Y and forming the product 

[Z] = IY][X] _1 

The subroutine DEIGEO solves the matrix problem of 
equation (3-4) . Subroutine MSET is called to set up matrix 
X. CDKTIN is then called to invert [X]. CDHTIN was 
obtained by modifying the IBM library routine CMTRIN to make 
it applicable to double precision matrices. Next, BMSET is 
called to set up matrix Y using space-saving band storage. 
The two matrices are multiplied using subroutine MULDBM . 
Since all the programs available for solving the resulting 
eguation require that the real and imaginary part of the 
matrix be presented in separate matrices, subroutine DSPLIT 
is called to accomplish this. The programs used to solve 
for the eigenvalues are EHESSC and ELRH1C which are 
available through the International Mathematical and 
Statistical library. These subroutines reduce the matrix 



(3-8) 



(3-9) 
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into the complex upper Hessenberg form and then solve for 
the eigenvalues. Another program which can be used to solve 
for the eigenvalues is EISPAC (Eigensystem Subroutine 
Package) which is available from the IBM library. This 
program was also used to solve this problem and gave results 
identical to the results obtained using the I.M.S.L. 
routines. DEIGEO solves the problem twice, once for the 
even eigenvectors and once for the odd eigenvectors. The 
results are then passed back to the calling program. 

Three main pregrams have been written. Program #1 
takes values of a and /3 as input, calls DEIGEO, and outputs 
a list of the eigenvalues, along with the stability for each 
eigenvalue, determined using equation (2-35). A printer 
plot is also output. Program #1A is the same as program #1 
except that a plot on the Calccrap plotter is output using 
subroutine DRAW. This output is used in this paper to 
present the graphs of the eigenvalues. Program #2 was used 
to determine the shape of the stability curves. This 
program performs five functions depending on the value of 
the ccntrcl parameter MODE. The stability curves determined 
are all two-dimensional plots of <X (the imaginary part of 

a) versus the Reynolds number, holding /B , /5 , a , the 

I R R 

velocity, and the stability constant. 



For MODE= 1 , the stability for a given set of input 
parameters, a, /3, Reynolds number (REY) and velocity (VEL) 
is determined and printed. Equation (2-35) is used to 
determine the stability. 

For M0DE=2 a point on the line of constant stability 
with an input value of STAB is determined for a given input 
value of REY. Referring to Figure 3-4a, the initial guess 
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has a value of and KEY, and corresponds to point 1. A 

second point (2) is determined by adding DAI to the initial 

value of a_^. The stability at both of these points is 

determined by calling the subroutine EIGLS which determines 

the value of the stability of the least stable eigenvalue 
for a given set of input parameters, a, /? , REY, and VEL. A 



linear 


interpolation 


on 


a is made 

I 


to determine a third 


guess 


(3) for which 


the 


stability 


is also 


calculated. 


Subroutine DFIT2 is 


then 


called to 


make a 


secon d-order 



approximation of the value of a_^ which corresponds to the 
desired value of the stability (STAB) . This final value is 
printed . 

For MGDE=3, the same operations are performed as for 
M0DE=2, except that REY is varied instead of a . Figure 

3-4b shows this. The initial guess is point 1, DREY is used 

to find point 2, a linear interpolation to find 3 and a 
second-order fit to determine a final estimate, which is 
printed. 

For M0DE=4, the point of minimum REY on the neutral 

stability curve is determined. Two points (1b and 1c) are 

determined from the initial guess (la) , by adding and 

subracting DAI frcm a . Their stabilities are determined by 

I 

calling EIGLS and a second order fit determines the least 
stable point (1) for the initial value of REY. A new value 
of REY is picked by adding or subracting DREY as 
appropriate. This results in point 2a. The points 2b and 
2c are found by adding ±DAI. The stability of these three 
points is determined, resulting in the least stable point 
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(2) for this value of REY. A third point (3a) is determined 
using a linear interpolation from points 1 and 2, two points 
(3b and 3c) are chosen by adding ±DAI and the stabilities of 
these points results in the least stable point (3) . 
Finally, a second-order fit is made using points 1, 2, and, 
3 to find the point corresponding to zero stability. This 
point will be the point of minimum Reynolds number on the 
neutral stability curve. 

For M0DE=5, the point of minimum stability is determined 

for a set of input values a , /? , 3 , and VEL. The 

R I R 

procedure is similar to that for H0DE=4. An initial guess 

(la) and two points (1b and 1c) are used to find point 1. 
This procedure is repeated twice at values of REY determined 
by adding or subtracting DREY from the initial guess. 
Points 1,2, and 3 are then used to approximate the point of 
minimum stability using subroutine DFIT2 . 

B. CYLINDRICAL COORDINATES 

Programs to solve the two cases, n=0 and n=1, were 
written separately. The method of finite differences is 
used in each case. All subroutines used in the cartesian 
case are used in the cylindrical case with the exceptions of 
the functions returning the coefficients and the subroutine 
setting up the matrices. 

Functions CHM1E1, CHM2E1, CHM1E2, and CHM2E2 return the 
values of the coefficients of the velocity vector potential 
component h for the five adjacent mesh points in the central 
difference representation of the equations. This is done in 
a manner exactly the same as the functions CHM1E1 and CHM2E1 
in the cartesian case. The ,, M2" means those terms which 
also contain the eigenvalue y, and are used to set up the 
matrix Y. The "HI" means those remaining terms which do not 
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contain y, and are used to set up the matrix X. "El" means 

that the terms come from the linear combination of T arid F 

x e 

given in equation (2-14) and which are written out as the 

top row of the matrices of (2-15) through (2-22) . "E2" 

means that the terms come from the equation T , and which 

r 

are the bottom row of the matrices of (2-15) through (2-22) . 



Similarly, functions CGM1E1, CGM2E1, CGM1E2, and CGM2E2 
return all the values of the finite difference coefficients 
of the component g (r) , and CFM1E1, CFM2E1, CFM1E2, and 
CFM2E2 return the values of the coefficients of f (r) . 

Subroutine I1SET2 performs exactly the same function as 
MSET, with the exception that MSET2 allows the matrix to be 
positioned with a starting value ether than (1,1). This is 
necessary because for cases other than n=0, two equations 
must be solved simultaneously and matrices must be placed 
adjacent to each other to accomplish this. Also, MSET 2 
forms a mesh which includes the point r=0, in contrast to 
MSET which had the center of the channel fall between mesh 
points. The boundary conditions must be written 
appropriately. MSET2 sets the correct boundary conditions 
at the v/all r=1, but does not set any boundary conditions at 
r=0 since they are different for each case. Thus, the 
boundary conditions at r=0 must be set by the calling 
program. 

1 . Axisymraet ric liPH ( n~ 0 ) 



Program #R1A solves the first equation for the 
function h (r) . As shown in part II, section C, the boundary 
conditions can be written 



h(1) = 0 



(3-10) 



Dh ( 1 ) = 0 



(3-11) 
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h (0) = 0 



(3-12) 



D2h (0) = 0 (3-13) 

All even derivatives at the origin are zero, but only (3-12) 
and (3-13) are needed to satisfy the difference equations at 
r=0. The above equations imply that the boundary mesh 
points at r=0 and r=1 are left out of the difference 
equations. When the virtual point outside the boundary r=1 
appears, its coefficient is added to the first mesh point 
inside the boundary. When the virtual point outside the 
boundary r=0 appears, the negative of its coefficient is 
added to the first mesh point inside the boundary. The 
boundary conditions at r=1 are set by subroutine MSET2 while 
the boundary conditions at r=0 are set by the main program 
#R1A. The rest of the solution proceeds exactly the same as 
the method used for the cartesian case discussed in section 
A of part III. 



Program #R1B solves the second equation for the 
function g (r) . The boundary conditions, from part II 
section C, are 



g(1) = 0 



(3-14) 



g (0) = 0 (3-15) 

In addition to (3-15), all odd derivatives of g are zero at 
r=0, but cnly (3-15) is required to satisfy the difference 
equations at r=0 since the highest derivative of g is the 
second derivative and the result is a tri-diagonal matrix. 
Again, the solution procedure is identical to the cartesian 
case. 



Both programs #R1A and #R1B output a list of the 
eigenvalues and the stability for each eigenvalue. Printer 
plots of the eigenvalues are also output. 
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Figure 3-5 Finite Difference Matrix for Pipe Flow 



2 . 



njJ 

For this case the two equations must be solved 



simultaneously. This is done by writing two sets of n 
equations for the set of 2n vciriables including n unknown 
values of f and n unknown values of g. The function h is 
set identically to zero as described in section II, part C. 



The boundary conditions are 
f(1) = 0 


(3-16) 


Df (1) =0 


(3-17) 


9(1) = 0 


(3-18) 


f (0) = 0 


(3-19) 


D2f (0) = 0 


(3-20) 


9(0) = 0 


(3-21) 


The matrices are set up as shown in figure 


3-5 using 



subroutine MSEI2 four times. The upper left section is set 
up using function CFM1E1. CGM1E1 is used for the upper 
right section, CFM1E2 for the lower left and CGM1E2 for the 
lower right. A second matrix is set up using an identical 
procedure, but for the ,, M2 U coefficients, that is, the 
coefficients which contain y. The procedure is then to 
invert this second matrix, multiply it with the first and 
solve for the eigenvalues. Program #R2 carries out the 
solution and then lists the eigenvalues found, with their 
stability. A plot is also output. 
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IV. RESULTS 



A. COMPARISON WITH PREVIOUS RESULTS 

1 . Neutral Stability Curve 

The neutral stability curve was calulated for the 

two-dimensional, sinusoidal case, that is, for a = )9=/3=0. 

R I R 

Figure 4-1 compares the results with those of Grosch and 

Salwen (1968) . Using n=6 0 (60 mesh points across the half 
channel) gives almost the same accuracy as the published 
results. The computation time for finding the eigenvalues 
using n=60 was over one minute (on an I.B.H. 360), compared 
to a computation time of approximately 10 seconds for n=30. 
For this reason a mesh size of 30 was used to obtain all 
results presented in this paper except where noted. 

2. Squire* s Theorem 

Squire* s theorem can be verified numerically. As 
expected, the rotation of three-dimensional sinusoidal 
disturbances onto the two-dimensional plane always has a 
stabilizing effect. From Sguire's theorem, it follows that 
the addition of a transverse sinusoidal component to a 
two-dimensicnal sinusoidal wave is absolutely stabilizing in 
that the minimum Reynolds number on the neutral stability 
curve must increase for increasing /3^ . Figure 4-12, which 

plots the minimum Reynolds number for various values of 

and B , shows this to be so. For a =0, the neutral 
I R 

stability curve does not seem to have another lobe (labelled 

e-f-g in figure 4-10a), as the curves do for a *0. For the 

R 

cases for which Squire's theorem applies, a =0, figure 4-12 
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Squire's theorem 



shows that increasing /3 is stabilizing. 

is not applicable to the other cases, where a *0. 

R 



B. STABILITY CURVES 

The stability criterion used is discussed in section 

II-D. It is based on how an observer in a moving coordinate 

system would experience the combination of exponential 

growth in space and in time of the perturbations. The 

stability is given by /' in equation (2-28) . 

R 

, yt = y + a U (2-28) 

R R R 

The stability observed from a coordinate system moving with 

velocity U is said to be stable, neutral or unstable 

according to whether v' is less than, equal to, or greater 

R 

than zero. The velocities of interest are those of the 
fluid particles, which vary from 0 to 1.5. 

The neutral stability curves for different values of a , 

R 

with /3 =0, are shown in figure 4-2. The stability criterion 

used is for a fixed coordinate system where U=0 in equation 

(2-28) . This is the same as the instability observed from a 

ncnmoving frame of reference. The structure of the curves 

for a ^0 appears to be qualitatively different from the 
R 

curve for a =0 in that (referring to figure 4-10a) the lobe 
R 

e-f-g is net present on the neutral curves for a =0. 

R 

Figure 4-3 shows the effect on the neutral curves of the 
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addition of a 0 component. It appears that only the upper 

lobe (labelled b-c-d in figure 4- 10a) of the curves is 

affected. The stability is y* from equation (2-28) , for a 

R 

fixed coordinate system (U=0) . 



Cases involving nonzero values of ft were not computed. 

R 

This was primarily because of a lack of time, but also 
because a hypothetical case of this kind might well be 
difficult if not impossible to realize experimentally. 

Curves of constant stability for a =0 and a =-.04 are 

R R 

shown in figures 4-4 and 4-5. It is important to note that 

equation (2-28) can be used (for nonzero values of a ) to 

R 

transform values of calculated stability as seen by one 

reference frame, into different values of stability as seen 

by another reference frame moving with a different U. 

Conversely, if a curve of constant stability is given for one 

reference frame, it is possible to find another reference 

frame which sees that curve as having any desired stability, 

such as neutral stability, for example. The velocity of 

this frame cf reference is meaningful only if it falls 

within the limits of the fluid particle velocities, 0 to 

1.5. Thus, for a =-.04, the neutral curve for a coordinate 

R 

frame moving with the mean velocity, U=1.0, corresponds to 

the curve of constant stability for y =.04 as seen from a 

R 

fixed coordinate frame. Correspondingly, the neutral curve 

for the maximum velocity, U=1.5, corresponds to the curve of 

constant stability for y =.06 for U=0. 

R 
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c. 



EIGENVALUES 



The method of finite differences was used to solve the 
differential eguations (described in section III) . A mesh 
of n finite difference grid points was placed across the 
half channel. 

The effect of mesh size on the eigenvalues is 
demonstrated in figures 4-6a,b ,c, and d. From these graphs 
it can be seen that the effect of increasing the accuracy of 
the solution by increasing the number of mesh points is to 
stretch out the eigenvalues. For high enough n, the 
eigenvalues form a "Y" lying on its side with the branches 
near the y axis. The other two branches, which move away 

from the y^ axis with increasing n, seem to be "false" 

eigenvalues which can be made ever increasingly stable by 

improving the accuracy of the solution. A true, qualitative 
picture of the eigenvalues is as shown in figure 4- 10b, with 
the line of eigenvalues extending indefinitely to the left. 
Note also that in all the eigenvalue plots shown, the value 
of y corresponding to the base of the "Y" formed by the 

eigenvalues is egual to the value of a . 



The effect of Reynolds number is shown by figures 
4-7a,b,c,d. It is apparent that changing Reynolds number 
has an effect similar to that of changing n. This seems to 
indicate that, for a given mesh size, lower Reynolds numbers 
yield more accurate eigenvalues than higher ones. This is 
shown to be true later in this section. 

The effect of a on the eigenvalues is shown by the 

R 
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series of plots in figures 4-7b, 4-8a, and 4-8b. The effect 

of decreasing a is to make the flow less stable. Note that 

R 

both branches of the eigenvalue plots become unstable. No 

plots of positive a are shown, but the results are to shift 

R 

the eigenvalues to decreasing y , that is, stabilizing. 

R 

The effect of is shown in figures 4-9a, 4-7b, and 

4-9b. Increasing a tends to increase the f regu ency { y ) , 

I I 

that is, shorter wavelengths are associated with higher 
frequencies . 

The significance of the different parts of the stability 

curves is shown in the series of figures, 4-11. The neutral 

stability curve for a =-.04 is shown in figure 4-10a with 

R 

the parts of the curve labelled. Figure 4- 10b shows a 

typical eigenvalue plot for large n with the branches 
labelled. Eigenvalues for the points labelled A, B, C, D 
and E, are plotted in figures 4-1 1a,b,c,d, e. From these 
plots it is observed that the traditional neutral curve, 
represented by the instability of point B, is associated 
with branch LF of the eigenvalues corresponding to lower 
frequencies. The curve a-b d-e that cuts vertically through 
the upper lcbe is associated with the instabilities of 
points A and D. These points have instabilities on the 
other branch (HF) of the eigenvalues, the one associated 
with higher frequencies. Figure 4-1 1c illustrates that at 
point C, both branches of the eigenvalues are unstable. At 
point E, figure 4-11e shows that the two branches, HF and 
LF, of the eigenvalues have almost disappeared, which seems 
to relate to the sudden decrease in stability of the lower 
lobe e-f. 
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D. 



STABILITY 



The instability of the top lobes (labelled b-c-d in 

figure 4-10a) of the neutral curves of figures 4-2 and 4-3 

is gualitatively different from instabilities elsewhere. 

For this reason, and because these lobes are directly 

analogous to the neutral curve for a =0, this region of the 

R 

curves is considered separately from the rest of the 

instabilities. Each of these lobes has a point of minimum 

Reynolds number. These are plotted in figure 4-12 for 

different values cf a and B . Once again, the stability 

R I 

criterion is based on equation (2-28) for a fixed coordinate 
system. 



Each of the upper lobes (b-rc-d) also has a maximum 
instability associated with it similar to the maximum 
instability shown in figures 4-4 and 4-5. These values are 
shown in figure 4-13 which plots the coordinates (a vs. Re) 

of the maximum instabilities. Only the case of B =0 is 

I 

considered. Note that the values of a plotted only go to 

R 

-.04. This is because the instabilities associated with the 



part of the neutral curve labelled a-b d-e in figure 4-10a 

overwhelm the instabilities of the lobe b-c-d as a gets 

R 



increasingly negative. 



If the values of the maximum instabilities are plotted 

vs. a , as in figure 4-14, they form the curve labelled 
R 
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U=0. The values c£ stability of figure 4-13 correspond to a 

ncnmoving coordinate system. Using equation (2-28), the 

maximum instabilities for nonzero velocities are plotted. 

The lines appear to be almost straight. Note that for a 

velocity cf approximately .41, the line of maximum 

instability is approximately horizontal. For velocities 

greater than .41 the line intersects the a axis at some 

R 

negative value cf a . For a below this value, the flow is 

B R 

stable. Assuming that these lines are straight, values of 

velocity less than .41 will make the flow unstable at 

Reynolds numbers as low as desired by choice of a suitably 

negative value of a . Although not yet verified 

R 

numerically, it seems possible that for velocities greater 

than .41 there may be a critical value of the Reynolds number 

below which no instabilities are found. Curves of constant 

stability for positive a must be plotted to see if this is 

R 

so. Cases of (3 *0 were not considered due to lack of time. 

I 



The most critical portion of the neutral curve occurs at 
values of approaching zero, that is, for long 

wavelengths. The stability of points along a line close to 

the Re axis (at a value of a =.01) was determined for values 

of a ranging from -.02 to -10.0. Equation (2-28) was used 
R 

to translate these values of the stability found for a 

ncnmoving coordinate frame of reference to the velocity of a 
coordinate axis for which the stability was exactly neutral. 
The resulting curves are plotted in figure 4-15. The region 
of instability lies below and to the right of the curves. 
It can be seen from the graph that there is no minimum 



53 



Reynolds number associated with these instabilities. In 

particular, a value of a =-10.0 seems to be completely 

R 

unstable for all velocities and Reynolds numbers. 



E. PIPE FLOW 

The results for the numerical solution of the pipe flow 
eguations seem to indicate that there is an error somewhere 
in either the analysis or the numerical solution technique. 
All three programs, R1A, RIB and R2, produce eigenvalues 
which become increasingly unstable for decreasing Reynolds 
number. Further study is necessary to resolve this. 



F. NUMERICAL ACCURACY 

Figure 4-16 shows the neutral stability curves for two 

values of a for a mesh with 30 points. Portions of the 
R 

curve with n=60 are also plotted. These graphs support the 

observation made earlier that decreasing Reynolds number 
seems to be associated with increasingly accurate solutions. 

Calculations were done on an I.B.M. 360 computer in both 
single and double precision, that is, with computer word 
lengths cf 32 and 64 bits respectively. Some calculations 
were also dene on an X. D.S. 9300 with a 48 bit word length. 
With n=60, the eigenvalues obtained on the 9300 and in 
double precision cn the 360 agreed to four to six places. 
Those obtained in single precision on the 360 sometimes 
agreed to only the first place with the double precision 

values. Thus it appears that the 64 or 48 bit word length 

is sufficient for matrices of size n=60 while in single 

precision cn the I.B.M. 360, some eigenvalues have 

significant errors. 
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Grosch and Salwen (1968) 
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Figure 4-1 Comparison with previously published results. 
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Figure 4-2 Neutral stability curves for various 
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Figure 4-3 Neutral stability curves for various /3j 
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Figure 4-4 Curves of constant stability for 
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Figure 4-6a n=30 ^-1.6 Figure 
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Figure 4 - 6 c n = 50 ^ - 1 .6 Figure 
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Figure4-7a Re s 4000 J -1.6 Figure4-7b Re : 8000 
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Figure 4-7c Re = 12,000 
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Figure 4-8a a R =-.04 J -l.6 Figure 4-8b a =-.08 
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Figure 4-9a = .4 Figure 4-9 b 
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Figure 4-IOa Neutral stability curve for a R = -.04, plot for large 
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Figure 4-lla Point A Figure 4-llb Point 
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Figure 4- lie Point C J " 1 - 6 Figure 4-tld Point 
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Figure 4 - I le Point 
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Figure 4-12 Points of minimum Re on 
neutral stability curves 
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Figure 4-14 Values of maximum instability of the 

upper lobe vs. Oi for various velocities 
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Figure 4-16 Numerical accuracy of two stability curves 



V. 



CONCLUSIONS AND RECOMMENDATIONS 



The solution of the linearized vorticity transport 
eguation was sinplified considerably by the introduction of 
the velocity vector potential. This method has not been 
mentioned in previous publications, at least not in those 
with which the author is familiar. 



The effect of the number of finite difference mesh 
points on the eigenvalues is important for understanding the 
effects of the ether parameters. This is explained in the 
results section, and is believed to be information not 
previously published. 



The results presented here show that the critical 

Reynolds number for plane Poiseuille flow can be lowered as 

much as desired by proper selection of ether parameters. In 

particular, it is highly significant that a negative value 

of a , that is, streamwise decay in space, produces highly 
R 

unstable rates of growth in time. The lowering of the 

critical Reynolds number in this way provides a new basis 
for improving the agreement between theory and experiment. 



Two separate modes of instability were found 

corresponding to the two branches of the eigenvalue plots. 

The classical case of instability found for a =0 corresponds 

R 

to only one of these branches becoming unstable. It is 

significant that for each different value of a , this mode 

R 

has a maximum instability associated with it. For a fluid 
particle velocity of approximately .41, this maximum 
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An additional 



instability seems to be invariant with a . 

R 

mode of instability, corresponding to longer wavelengths, 

appears for nonzero a . No maximum value of instability is 

R 

present as in the first case. This mode has been shown to 

be unstable at all values of Reynolds number for 

sufficiently negative values of a . 

S 



The stability of a fluid particle has been shown to 

depend upon its velocity. The fact that negative a yields 

R 

the greater instabilities indicates that fluid particles 

\ 

with the lowest velocities, that is, those nearest the 
walls, will be the most unstable. This seems to agree with 
experiment. For example, from Schlichting (1968), the 

transition is "characterized by an amplification of the 
initial disturbances and by the appearance of 

self-sustaining turbulent flashes which emanate from fluid 
layers near the wall along the tube". Further comparison 
with experiment is recommended. 



The eigenvectors of the solutions to the vorticity 
transport eguaticn were not investigated. Those 
corresponding to the least stable eigenvalues should be 
studied. 



Since satisfactory results were not obtained for pipe 
flow, apparently because of an error, an independent study 
of this case should be made, including the boundary 
conditions and the numerical methods. It seems likely that 
the solutions will be qualitatively similar to those 
obtained for plane flow. If this be true, then 
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instabilities should be found for negative values of a . 

R 

This would be significant because it is generally agreed 

that for a =0, pipe flow is stable to infinitesimal 
R 

perturbations. 
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APPENDIX A 



DERIVATION OF 

VORTICITY TRANSPORT EQUATION 

The flow of an incompressible fluid with constant 
viscosity is governed by the Na vier-Stokes equation 

w 2 v - (V • v) v - vp / p - av /at = o (A-i) 

d d d d d d 

where the subscript d indicates the dimensional quantities 

of the velocity vector, pressure, and time (V ,p , and t ). 

d d d 

The kinematic viscosity 1/ and the density p are constants. 

This equation is nondimensionalized using a reference 
length L equal to the half width of the channel in plane 
flow or the radius of the pipe in pipe flow. The reference 
velocity is the average velocity of the laminar flow V 

avg 

The resulting eguation is 

i_v 2 v _ (V.v) V - vp - av/at = o (a-2) 

■Re 

where Re f IV / v is the Reynolds number. This flow is 
avg 

also governed by the continuity eguation. 



v*V = 0 . (A-3) 

Other than equation (A-1), all the equations presented in 
this paper are in ncndimensional form. 

The vorticity transport equation is obtained by taking 
the curl (vx) of the Na vier-Stokes equation. The definition 
of the vorticity is 
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fl = VX V 



(A-4) 



and the following vector identities are applicable; 



v*vp f 0 



(A-5a) 



VX[(V«v)V] = VXt V (V«V) /2-VX (VXV) ] = -vx(Vxft) (A-5b) 



vx(v 2 V) = v z £l 



(A-5c) 



vx(3V/<?t) = a&/3t 



(A-5d) 



Combining equations (A-4 ) , (A— 5) , and the curl of (A-2) , the 

result is 



Next, it is assumed that the flow is made up of the 
steady-state, laminar solution for a given flow geometry and 
a small perturbation flow which is added to the laminar 

flow. The velocity and vorticity become, respectively, V+v 

and il+a) where v and u> are the perturbation guantities. It 

is easy to show that the perturbation flow must satisfy 
continuity independently. 

W = 0 ( A— V) 

It can also be shewn that 



Substituting V + v and il+w into the vorticity transport 



i v 2 X2 + vx (vx&) - ail/at = o . 

Fe 



(A-6) 



w - vx V 



(A- 8) 



78 



equation and subtracting the equation for the unperturbed 
flow, the result is 



1_V 2 o» + 7X( VXtu-^XV+vxw) - = 0 . 

He "31 



(A-9) 



This eguation is linear in the perturbation quantities 

except for the second order term v*o>. if the perturbation 

is small, this term can be neglected, linearizing the 
equation. The result is the linearized perturbation 

vorticity transport eguation. 

1 v 2 <*> + vx(Vxw-,Qxv) - =0 (A- 10) 

' He at 

A second form of this equation is derived by expanding 
the curl cf the cross product terms and noting that v*V = 
v • v = = vco - 0. The result is 

(1/Re) v 2 ^ + (o»#v) V - (V»v)cu + (v«v)fl 



A third form of eguation (A-10) is obtained using the vector 
identity 



(,fr*V) V - doJ/dt = 0 . 



( A- 1 1 ) 



v 2 a> = v(v»w) - vx(vxoj) 



(A- 12) 



and noting that a> is non-divergent 



-1 VXVXfl^ + VXVXaJ - VX^XV - 3a» = 0 
He St 



(A- 13) 
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APPENDIX B 



DEPENDENCE BETWEEN 
COMPONENTS OF THE 

PERTURBATION VORTICITY TRANSPORT EQUATION 



The vorticity transport equation 

-1_vxvxa> + VXYxw - vx&xv - 9 oj = 0 ( A— 13) 

Be 

appears superficially to represent a set of three, coupled, 
simultaneous equations. This equation actually has only two 
independent conditions. To show this, equation (A-13) is 
first simplified by the introduction of the velocity vector 

potential W defined by 

v = vxW . (2-7) 

In Appendix D it is shown that this results in a vector 
equation which contains only the three components of W as 
its unknowns. 

-1 vxvxvxvxW + vxVxvxvxW - vxflxvxW 
Be 

- vxvx9w = 0 { D — 1 ) 

Furthermore, the following form of solution is assumed. 

X 

W(x,y,z,t) = £ if (y) + jg (y) + kh (y) ]e (2-9a) 



where 



X = ax + /3z + y t • (2-10) 

Consider equations (D-1) expressed in abbreviated form by 
the statement 
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r = < r > = 0 
y 




( B 1 ) 



which may be expressed as three separate equations. 




(B-2a) 




(B-2b) 




(B-2c) 



Equations (B-2) are written for cartesian coordinates. The 
components T , T and T happen to be the three components 



of a vector obtained by applying the curl operator to the 

momentum transport equation. It is a vector identity that 
the curl cf any vector is nondivergent and thus 



With the assumed form of W given by equations (2-9a) and 
(2-10) and the fact that all these equations are linear in 
W, (B-3) becomes 



Equation ( B— 4) has been verified for the resulting form of 

the eguation developed in Appendix D (equations ( D— 46) 
through (D-55) ) by differentiating the second eguation , 
multiplying the first by a and the third by and adding 

the three equations thus obtained. All terras are found to 
cancel out. From eguation (B-4) it may be inf erred that it 



x y 



z 



^_(T) + 3 (D + i_(T) = 0 

dx x 3y y 3 z z 



(B-3) 




(B-4) 
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is sufficient tc satisfy the second of equations (B-2) and 
any linear combination of the first and third. If the 
second is satisfied, then 

T = 0 (B-5) 

y 

and it must be true that 

|_ ( T ) = 0 . (B-6) 

y 

Eguaticn (B-4) reduces to 

a r + py = 0 . ( B— 7 ) 

x z 

As discussed in section 2, the proper choice of a second 
condition is 



-pr + oT = 0 

X 2 



( 2 - 11 ) 



Equations (B-7) and (2-11) may be expressed in matrix form 
as 



a P 
‘P a 




(B-8) 



If the determinant is non-vanishing, that is, if 

a 2 + P z * 0 , (B-9) 

then the only solution of equation (B-8) is 



(B-10) 



Thus, equations (B-5) and (B-10) verify that equation (B-1) 
will be satisfied by the two equations (B-5) and (2-11). 
That is, these two conditions are sufficient. 



For pipe flow in cylindrical coordinates, equation (D-1) 
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may be expressed similarly to (B-1) as 



r = 



s r 



> = o 



e J 



( B- 1 1 ) 



which are three separate equations. 



r = o 

x 



(B-1 2a) 



r = o 

r 



(B- 12b) 



r Q = o (B- 12c) 

As in the cartesian case, equation (B-11) must be 
nondivergent. 

) + Ji_ (r T ) ♦ li ( T ) = 0 (B-1 3) 

ax x rar r raS 0 



With the assumed form of W given in equations (E— 1) and 
(E-2) , eguaticn (B-13) becomes 



or + jd (r r ) + £3 r = o 

x r r rQ 



(B-14) 



Eguation (E-14) has been verified by performing the 

indicated operations on T , T , and T . All terras cancel 

x r e 

out. 



From eguation (E-14) it may be inferred that it is 

sufficient to satisfy equation (B-12b) and any linear 

combination of T and T . If the second is satisfied then 

x G 



r = o . 
y 



(B-15) 
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Therefore 



ID (r D = 0 (B- 16) 

r r 

for any value of r except at r=0, where it is true in the 
limit as r approaches zero. Equation (B-14) then reduces to 

aT '+ §T =0 (B-17) 

x r 8 

As discussed in section 2, the proper choice of a second 
condition is 

-§ T + aT =0 (2-14) 

r x 9 

These two equations may be expressed in matrix form as 

’ a £ 
r 

-2 a 
r 

and, as before, equation (B-11) will be satisfied under the 
condition that 

a z + * 0 . (B- 19) 

r " 2 
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APPENDIX C 



REDUNDANCY OF THE THREE COMPONENTS 
OF THE VELOCITY VECTOR POTENTIAL 
AND THEIR RELATIONSHIP TO THE STREAM FU NOTION 

1 . Cartesian Coordinates 

To find the relationship between the components of 
the velocity vector potential in cartesian coordinates, 

consider the definition of W. 





i j 


k 




i 


3 


k 






V = VXW = 


3 <9 

3x • ay 


d 

3z 


X 

e = 


a 


D 




X 

e 






f (y) g{y) 


h (y) 




f 


g 


h 




(C-1) 


which has the components 
















u (y) = Dh - 


■ $ g 














(C-2a) 


v(y) = £f - 


■ ah 














(C-2b) 


w(y) = ag - 


• Df . 














(C-2c) 


The components 


of the 


velocity. 


u. 


V 


/ 


and 


v, are not 



independent but are related by the continuity equation, that 
is, by the condition of non-divergence. 



v» 



u = au + Dv + ySw = 0 



(C-3) 



Let f (y) be identically zero. Then 
u(y) = Dh - /3g 

v (y) = -ah 



(C-4a) 

(C-hb) 
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v(y) = ag . 



(C-4c) 



Thus the components of the velocity may be represented by 
the two components, h and g, of the velocity vector 
potential without any loss of generality other than that 
imposed by the condition of non-divergence. If g (y) is set 
to zero, the following set of equations result. 



u(y) = Dh (C-5a) 

v(y) = 0f - ah (C-5b) 

w(y) = -Df (C-5c) 

And if h (y) is identically zero 

u (y) = -£g (c-6a) 

v(y) = (c-6b) 

w(y) = ag - Df . (C-6c) 



In the classical case of plane perturbations the 
parameter /? is zero. Then the 3-dimensional equations, 
(C-2) , reduce to 

u (y) = Dh (C-7a) 

v (y) = -ah (C-7b) 

w (y) = ag - Df . (C-7c) 

From equations (C-7) it can be seen that, for the case of 

/3=0, letting h(y)=0 cannot be allowed because both u(y) and 
v (y) are forced to be zero. Since it is desirable to have 
the three-dimensional representation of the problem 
reducible to the classical two-dimensional case, it is 
stipulated that the function to be eliminated not be h(y) . 
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Note that the classical case of two-dimensional 
plane perturbations corresponds, in the present notation, to 
setting f (y) and g (y) egual to zero and retaining h (y) 
alone. This component of the velocity vector potential 
corresponds to the classical definition of the stream 
function for two-dimensional flow. From equations (C— 7 ) 

u(y) = Dh (C-8a) 



v (y) - -ah (C-8b) 

w(y) =0 (C-8c) 

for the classical two-dimensional flow. 

2. Cylindrical Coordinates 



The relationship between velocity and the vector 
potential in cylindrical coordinates is similar to that 
relationship just developed for cartesion coordinates. The 

definition of W in cylindrical coordinates is 



v = vxw 



Je _1e 
r x r r 



a D 

f g 



e 

e 




which has the components 



u(r) = 1D(rh) - £g 
r r 



(C-9) 



(C-lOa) 



v (r) = Qi - ah (C-lOb) 

r 

w (r) = a g - Df (C-IOc) 

As before, the components of velocity are related by the 
continuity equation 
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(C-11) 



v*u = a U + _1D (rv) + = 0 . 

r r 



It is easily shown that the velocity components in 
cylindrical coordinates also can be represented by two 
components of the velocity vector potential, without loss of 
generality other than that imposed by the condition of 



nondivergence. In the case of 
eguations become 


P= 0 and h (r) = 0, the 


u (r) =0 


(C-12a) 


v (r) =0 


(C-12b) 


w (r) = ag -Df 


(C- 12c) 



which cannot be allowed. Boundary conditions in cylindrical 
coordinates impose the following restriction on P. 

P = in, n = 0,1,2,... (C-13) 

The case of n=0 is considered separately from n#0. For n=0 
it is necessary to retain h (r) and set either f (r) or g (r) 
egual to zero. For the other cases, n#0, h (r) may be the 
component of the vector potential which is set to zero. 



Note that n=0 corresponds 

perturbations. For f (r) =0 


to axially symmetric 


u (r) = ID (rh) - ^g 
r r 


(C- 14a) 


v (r) = -ah 


(C - 14b) 


w (r) = ag . 

For g (r) = 0 


(C- 1 4c) 


u(r) = iD(rh) 
r 


(C-1 5a) 
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(C- 1 5b) 



v (r) = &£ - a h 
r 

w (r) = -Df . 

And for h (r) =0 

u (r) = -gq 
r 

v (r) = £f 
r 

w (r) = a g -Df . 



(C- 1 5c) 

(C- 16a) 
(C- 1 6b) 
(C- 16c) 
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APPENDIX D 



DERIVATION CF COEFFICIENTS FOR 
PLANE POISEUILLE FLOW 

The hasic eguations are 



-1 VXVXW + vXVXu - vxllxv - du)/dt = 0 

I?e 


(A-13) 


OJ = VX V 


( A— 8) 


V = i3(1-y2) = iU (y) 


(2-3) 


= k3y 


(2-4) 


v = vxW 


(2-7) 



The vector functions W, v and a> are assumed to have 



solutions of the form of eguations (2-9) . 

X 

w (x,y, z,t)= [ if (y) + jg (y) + kh (y) ]e 


(2-9a) 


- - - X 

v(x / y / z / t) = [iu (y)+jv (y) +kw (y) ]e 


(2-9b) 


- - - X 

w(x / y,z / t) = £ i f(y) + ji?(y) +k C(y) ]e 

where 


(2-9c) 


X ~ ax + /8z + yt 


(2-10) 



Equation (A-13) can be written in terms of only W by 
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substituting for oj using equation (A- 8) and then 
substituting for v using equation (2-7) . The result is 



-1_vxvxvxvxW + vxVxvxvxW - vxfUvxW 
He 



- VXVXdW = 0 
HI 



(D-1) 



From the assumed form of the solution for the velocity 

vector potential W, given in equations (2-9a) and (2-10), it 
can be seen that the partial derivatives 5 /3x, d /Hz, and 
5/<?t become multiplication by a , J3 and y, respectively. The 
partial derivatives 5/ay r dz/Hy , etc. will be represented by 
the operator symbols D, D 2 , etc. 



The easiest method for determining the form of the 
coefficients of eguation (D-1) is to represent the operators 

as matrices multiplying the vector W and its derivatives. 

Brackets are used to indicate a 3 by 3 matrix. Using braces 



to enclose the components 
definitions are assumed; 



f (Y 



w = 



a vector, the following 



(D-2) 




DW = 



Higher order derivatives are defined similarly, 



(D-3) 



It is necessary to determine the matrix operator [CP] 
which can be used to replace the cross product operator, 
that is, a matrix [CP] which satisfies the relation 



[ CP ]W = tx[s]W 



(D-4) 
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where t is an arbitrary vector and [s] is an arbitrary 
matrix. If the components of t are t , t and t and the 

12 3 



components of [s] are s , then 

ij 



ft 1 


— 


1 




k 




t 




\ 3j 





s 


S 


S 


1 1 


1 2 


l 3 


s 


S 


S 


2 1 


2 2 


2 3 


S 


S 


S 


3 l 


32 


3 3 



'f' 



9 

h 



( D— 5 ) 



The matrix [CP] may be obtained by multiplying vector W by 
matrix [s], then finding the cross product of the resulting 
vector with vector t and, finally, collecting terms in f, g, 
and h to separate out the vector W. The result is 



[CP] = 



(t s -t s ) (t s -t s ) (t s -t s ) 

2 31 3 21 2 32 3 22 2 33 3 23 

(t S ~t S ) (t S ~t S ) (t S -t S ) 

3 11 1 31 3 12 1 32 3 13 1 33 

(t s -t s ) (t s -t s ) (t s -t s ) 

1 21 2 11 1 22 2 12 1 23 2 13 



(D-6) 



A matrix operator is also needed to replace the curl 
operator, vx. 

I j k 

3 3 3 

3x 3y 3z 

XXX 
fe ge he 

— X — X “ X 

= i (Eh-/%|) e + j(/3f-c(h)e + k(ag-Df)e ( D— 7) 
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This can he written 



VxW = [ A ]W + [ B ]DW 



(D-8) 



where [A] and [B] can be determined by collecting the proper 
terms. 



[A] 



0 -y8 0 

£ 0 -a 

_0 a 0 



(D-9) 



[B] 



0 0 1 

0 0 0 

-1 0 0 



(D- 10) 



This operation defined by equation (D-8) holds for the curl 

of any vector with the form (u (y) ,v (y) , w (y)^ e , that is, for 
any curl operation performed in equation (D-1) . 



To evaluate the third term in equation (D-1) take the 
cross product of Q, and vxW, that is, 
jjQx(£A]W + £ B jDW) . 

This can he evaluated using equations (D-6) and (D-5) . The 
result is 



ftxvxW = £C ]W + £ C ]DW 
1 2 



where 



[CJ = 



-3/Sy 0 3ay 

0 - 3/8y 0 

0 0 0 



(D-1 1) 



(D- 12) 
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[C 



2 



] = 



o 

0 



0 0 
0 3y 



0 0 



0 



Next operate with -vx. 



-vx ([C i ]W + [ C 2 ]) DW 

= - [A](rC i ]W+£C 2 ]DW) - [B]D([C i ]W+[C 2 ]DW) 
= - C A ]£C m JW - [A][C 2 JDW - [BjtCJDW 

- [BMcjw - [b][C 2 ]d 2 w - [B]D[C 2 ]DW 
= ( - £ A ][ C ^ ] - [B]D[C t ])W 

+ (-£Aj£C 2 3 - [BtfCJ - £B]D£C 2 ]DW 

- £E][C 2 ]D 2 w 



Note that D[ C ] means 9 [C ]. The coefficient of 

i 3y 1 



egual to the zero matrix, so (D-14) may be written 



-VXftxvW = fE^W + £E 2 ]DW . 

The matrices fE^] and £ E 2 ] may be found by simple 
multiplication and addition. 



[E 



l 



] = 



0 -3/32 y 0 

3/3 2 y 0 -3a/9y 

-3/8 3a/3y 3a 



(0-13) 



(D-14) 

D 2 W is 

(D-15) 

matrix 

(D-16) 
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[El = 



0 0 3(3y~ 

0 0 0 



-3/3y 0 0 



(D-17) 



To calculate the other terms, it is necessary to first 

find vxvxK 



VX (VX’W) = VX{[A]W + [BJDW) 

=[A][A]W + [ A ][ B ]DW + [ B ][ A ]DW 

+ [B]D[A]W + [B][B]D 2 W + [B]D[B]DW (D-18) 

This may te written 



VXVXH = [E^W + £E^]DW + [E s ]D 2 W 



(D- 19) 



where the matrices [E ], [E^] and [ ] are determined from 

the coefficient matrices of equation (D-18) . 



-£ 2 


0 


a/3 


0 


- (a 2 +/3 2 ) 


0 


r — 

p 


0 


-a 2 



[E 4 ] = 
[E s ] = 



0 a 0 

a 0 /3 

0/3 0 

- 1 0 0 

0 0 0 

0 0-1 



The term -vxvx^W can now be evaluted. 



(D-20) 



(D-21) 



(D-22) 
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( D— 23) 



-vxvx<5w = -yVXVXW 

3t 7 

= -y([E 3 ]W + [E^]DW + [E s ]D2W) 



To evaluate the term vxVxvxvxW, first find the matrix 
operator that is eguivalent to the cross product of V and 
vxvXW. Using eguations (D-19) and (D— 6 ) 



Vx (vxvxH) 



*<[E 3 



]W+[ E ]DH+[ E ]D2») 

4 5 



= £C 3 ]W + [C^]DW + [C ]D2W 



(D-24) 



where 



f c ,) 



0 

-a/3U 



0 

0 



0 

ct 2 U 



0 - (a 2 +/3 2 ) u 0 _ 



(D-25) 



ic.] - 



0 0 0 

0 -/3u 0 

_aU 0 /3U_ 



t c 5 ] - 



0 0 0 

0 0 u 



.0 0 0 . 

How operate with vx. 



( D— 26 ) 



(D-27) 
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VXVXVXVXW = VX (( C 3 ]W+[ ]DW4[ C g ]D 2 W) 

= (£A]+£B]D) ([c ]W+[C JDW +[ C ]D 2 W) 

^ ^ 5 

= £E 6 ]W + [E 7 ]DW + [E e JD2W + [E 9 ]D3H 



where 



tV 



'V 






a/3 2 U 3 (a 2 +/3 2 ) y -a 2 /3u 
0 a(a 2 +/3 2 ) 0 0 

- a 2 /3U 0 a 3 U 

-3ay -a 2 U -3/3y 

- a 2 U 0 -a/3U 

. 0 -a£u 0 

aU C O' 

OOO 

. 0 0 all. 



£E 9 ] = £0], the zero matrix. 



(D-28) 



(D-29) 



(D-30) 



CD-31) 

CD-32) 



To evaluate the term -1_vxvxvxvxW, first apply the curl 

Ee 

operator to vxvXW which is given by equations ( D— 19) through 
(D-22) . 

VX{ 7 XVXW) = ([ A ]+[ B ]D) ([ E 3 ]W+[ E^ ]DW+£ E g ]D 2 H) 

= £c ]w + £ C 7 3 DH + [ C 8 3 D2W + C C 9 3 D3W (d-33) 
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where 



- 



0 



P(a 2 *0 2 ) 

0 



“ fi(a 2 +/3 z ) 

0 - a(a 2 +/3 2 } 



0 

a(a 2 +£ 2 ) 

0 



[C 



7 



] = 



0 0 - (a 2 + /3 2 )" 

0 0 0 



(a 2 +/3 2 ) 0 0 



I°s’ - 



' 0 £ <f 

-0 0 a 



0 -a 0 



[C 



9 



3 = 



o o-i 

0 0 0 



1 



0 



0 



Finally, operate once more with the curl. 



vx ( vxvxvxW) 




+ [E 



= ([A]+[B]D) {[ C ]W+[C 1DH*[C ]D 2 W+[C 

O / o 

]W + [ E ]DW + [E i2 ]D 2 W + [E 13 ]D3W 
]D 4 W 

14 



where 



[E 




0 2 (a 2 +0 2 ) 0 - a0(a 2 *0 2 ) 

0 (a 2 +/3 2 ) 2 0 

- a0[a 2 +0 2 ) 0 a 2 (a 2 +0 2 ) 

» *4 



(D-34) 



(D-35) 



( D— 36) 



(D-37) 



] D 3 W ) 



( D— 38 ) 



(D-39) 
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[E ] = 
11 



0 - a(a 2 +/3 2 ) 0 

- a{a 2 +/3 2 ) 0 -/3{a 2 +/3 2 ) 



0 



-/3(a 2 +/3 2 ) 



[E ] = 
12 



a 2 + 2/3 2 o — aft 
0 {a 2 +/3 2 ) 0 



[E ] = 

13 J 



[E ] = 

14 



-a/3 

0 -a 0~ 
-a 0-/9 

0 -/9 0_ 

1 0 0 * 

0 0 0 

0 0 1 



2a 2 +/3 : 



The final eguaticn becomes 



-1 [E ]D^W - 1 [E ]D3?J + (-1 [E ]+[E ])D2^ 



He 



We 






DW 



+ (-1 £E in ]+[E ]+[E 1)W 

N Ee 1 0 6 1 ' 

- y([ E s ]D2W+[ E^ ]DW+[E 3 ]w) = 0 



which can be written 



(D-40) 



(D-4 1 ) 



(D-42) 



(D-43) 



(D-44) 



[M ]D*W + [M 3 JD 3 W + [ M ^ JD 2 W + [M^DW + [ M q ]W 

+ y([N 2 ]D2W + [ N ^ ]DW + [N^wJ =0 (D-45) 



where 
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-1 

Be 



[M 



4 



] = 



0 



0 

0 



0 

0 



0 



0 



-1 



Be 



[M 



3 



] = 



0 a 0 
Be 



a 0 B 
Be Be 




[« 2 ] = 



-£ 2 +T o a£ 

Be Be 

0 -1 {a 2 +£ 2 ) 0 

Be 

0 -a 2 +T 

Be Be 



[H^ = 
[H o ] = 

c» 2 ] = 



- 3 ay -aT 
-aT 0 

- 3 /Sy -/ST 

/S 2 T 
3 / 3 2 y 

a/ 0 T- 3/3 

1 0 o’ 

OOO 
0 0 1 



0 

-/ST 

0 

3a 2 y 
(a 2 +/S 2 ) T 

3a/3y 



-a/ ST 
-3a/Sy 

3a+a2T 



[N t ] = 



0 -a 0 
-a 0 -/S 



0 -/S 0 



(D-46) 



(D-47) 



{ D-48) 



(D- 49 ) 



(D- 50 ) 



(D-51 ) 



(D- 52 ) 
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( D— 5 3 ) 



£ 2 


0 


-a/3 


0 


{ a 2 +/9 2 ) 


0 


-a/3 


0 


a 2 



and T is defined by 

T = aU - 1_(a 2 +/3 2 ) = 3 a( 1-y 2 ) 
Ee 1 



1 

Ee 



(a 2 +/3 2 ) . 



(D-54) 



Finally, equation (D-45) 
yield the results below. 




can be divided through by e to 
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APPENDIX E 



DERIVATION CF COEFFICIENTS FOR 
PIPE FLOW 



The basic eguations are 

-1 vxvxvxvxw + vxVxvxvxW - rxilxvxw - vxvx aw =0 (D-i) 

le n 

V = e 2 ( 1-r 2 ) = e 0 (r) (2-5) 

xx 



£1 = e 4r (2-6) 

e 

and the velocity vector potential W has a form similar to 
that for the plane flow case. 

X 

W(x,r,0,t) = £ e f(r)+e g(r)+e h(r)]e (E— 1 ) 

x r 6 



where 



X = ax + @6 + yt . { E— 2) 

As in the plane flow case, the partial derivatives 3 /<3x, 
2/de and 3/dt become multiplication by a, $ and y, 
respectively. The partial derivatives <5/0r, d z /dz z , etc. 
will be represented by the operator symbols D, D 2 , etc. 



The nethod of solution for cylindrical coordinates 
proceeds exactly as that for the cartesian coordinate case 
presented in Appendix D. Definitions similar to { D— 2) and 
(D-3) are assumed. 



W 




X 

e 



(E-3) 
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DU 




X 



e 



(E-4) 



Higher order derivatives are similarly defined. The matrix 
[CP] defined by eguation (D-6) is valid for cylindrical 
coordinates and is used to evaluate cross products. The 
curl operation in cylindrical coordinates is found below. 




e 

e 



d 

33 

X 

rhe 



_ x — x — x 

= je (h+rDh-/3g^ e + Je ^Sf-rah^e + e^(ag-Df)e ( E— 5) 



This can be written 



vxW = [ A ]W + [ B ]DW 



(E-6) 



where 



£ A ] = 



-£ 

r 



i 

r 



[B] = 



0 a 

0 0 

0 0 

-1 0 



0 -a 



(E-7) 



(E-8) 



Following the method of Appendix D, the first step is to 
use eguaticn (D-6) to evaluate iixvxW. 



ibtvxii =-• JoJ * (£ A 3+c b 3DW) = [CJW + [C 2 ]d;; (E-9) 
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where 



C c t ) - 



- 4/3 0 4 ar 

0 - 4/3 4 



(E- 10 ) 



[c 2 ] = 



0 0 0 
0 0 4 r 

0 0 0 



t E— 11 ) 



Next, operate with -vx 



-VxiixvxW = ~([A]+[B]D) ([ C ]W+[ C 2 ]DW) 



= £E ]H + £E a ]DW + [ 0 ]D 2 W 



(E— 12) 



where 



£ E t ] 



0 -_ 4 / 3 £ 4/3 

r r 



4 / 3 £ 

"r 



0 -4 a /8 

4 a /3 0 



(E- 13 ) 



£E 2 ] = 



0 0 up 

0 0 0 

— 4/3 0 0 



(E- 14 ) 



Next, find rxvxW. 



VX(TXH) = (£ A )+£ E ]D) (£ A 3 W+£B 3 DW) 



= C e 3 3 h + £ E 4 3 DW + £E g ]D 2w 



(E- 15 ) 



where 
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-f! 




a 

r 




1 

r 


ce 3 ] = 


0 


- (a 2 +g 2 ) 
r 2 


r* 




r 

L_ 




~r^ 




-a 2 + 1 
F 2 




~-1 

r 


a 


o' 






[BJ = 


a 


0 


-£ 

r 








0 


§ 

r 


-1 

r 







C E ] 

5 



-1 0 0 

0 0 0 

0 0-1 



The term -vxvx^W can now be evaluated. 



-vxvx3w = -xvxvxH 
31 

= -x(ce 3 ]h + £E 4 ]DW + [E s ]D2W^ 



Next find VxvxvxH using eguations (E-15) and ( D— 6 ) 



Vx(vxvxH) = 




»([E ME lDWt[ E ]D 2 W) 

3 4 5 



where 



= C c 3 ] w + C c 4 3 dw + [ C 5 ] D2W 









-agu 

r 



& 

- (a 2 +£ 2 ) U 
r 7 



(a 2 -1 ) U 
r^ 

J3 U 
r* 



(E- 16) 



(E- 17) 



(E- 18) 



(E- 19) 



( E— 20) 



( E-21) 
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tc,] 



0 0 0 

0 -gU JU 
r r 



aU 



0 £U 
r 



(E-22) 



[c s ] = 



0 0 0 
0 0 0 

0 0 0 



(E-23) 



Now, operating with vx, the result is 



VXVXVXVXW = VX ([C 3 ]W + [C^ ]DW+[C 5 ]D2«) 



= £E 6 ]W + £ E ? ]DW + £E fi ]D 2 W + [E^D^W 



(E-24) 



where 



£ E ] 

o 



agl U 4r (a z +g z ) -a£0 -Hgg-ga z \l 
r 2 r 2 r r r~ 



-af£c 

r 



a(a 2 +^2) o 



ct^U 

r ? 



-aBU 
r 2 

(a3-a) u 
r 



(E-25) 



tV = 



aU-4ra -a 2 U -4/? 
r 



- a 2 0 
0 



0 -a£U 
r 



-a£u 

r 



aU 

r 



(E- 26) 



£E fl ] = 



aU 0 0 

0 0 0 

0 0 aU 



and the matrix £E ] is identically zero. 

9 



(E-27) 
For the final 



term, first calculate vxvxvxW. 
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VX(VXVXW) = ([A]+[B]D) ([E 3 JW+[E^ ]DW+[E s ]D2W) 



= [C fe ]W + £C ? ]DW + [C e ]D2W + [C 9 ]D3W 



where 



tV 



[ c 7 ] = 



0 


_4 + £(a 2 +£l) 

r 3 r r 2 


"l -J(a 2 +£ 2 ) 
r 3 r 


-£(a 2 +£|) 

r 


2 §f 


- a + a(a 2 +/3 2 ) 
r ? r* 


- 2 @ 2 
r 3 


a -a(az*||) 




2 |f 


0 




+ 1 
r* 






0 a 

r 






(a 2 +^|)-1 
£ r r r 


-a 0 

r 







tc e ] = 



s 

r 



-2 

r 



r 

1 

r 



-a 



0 a 
0 



£C 9 ] = 



0 0-1 

0 0 0 

1 0 0 



And, finally, 

VX(tfXVXVXii) = ([A] + [B]D) (£ ]W+[ C 7 ]DW + [C 0 ]D 2 H + [C ( 

= £ E ]W + £ E ]DW + [E ]D 2 W + [ E ]D 3 W 

10 11 1^ ^ ^ 



+ £ E ]D 4 W 

*■ 14 



(E-28) 



(E-29) 



(E-30) 



(E-31) 



(E-32) 



]D 3 W) 



(E-33) 



107 



Letting 



t = (a 2 +£2) 

r 2 



(E-34) 



the coefficients are 



[E ] = 
10 



2S|£ 

-a@t 

r 



- a -at 
r 7 r 

^2- a 2+t2 

r^ r 2 



-aB-aBt 

r 

-3 -2a2-3/3s+a2t 
r T i 2 r T 



(E-35) 



[E ] = 

li 



■3^ 2 +J +t 
r 7 r 7 r 

-at+ a 
r 2 

-aB 

r 2 



a -at 
r 2 

-^2+a£ 

r 7 r 



aB 

-£t+£ 
r r 7 



-3fi-fit 3 +a 2 -20 2 +t 
r 7 r r 7 r~ r 7 r 



(E-36) 



[*i 2 ] = 



-1 +t 
r 2 


-2a 

r 


-2§ 

r 


-a 

r 


t 


-2$ 

r 2 


-aB 

r 


2B 

r 2 


rf2-3 +t 

r 2 



(E-37) 



[E 3 

13 



2 

r 

-a 



-a 



0 - 



2 

r 



-fi 

r 

2 

r 



(E-38) 






1 0 0 
0 0 0 

0 0 1 



(E-39) 



The final equation/ with all the terms grouped together 

X 

has the same form as equation (D-43) and, dividing by e , 
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can be written in the form of equation (D-55) 







=0 



where 



£MJ = 



-1 0 
Be 



0 

0 



0 0-1 
Be 



£« 3 ] = 



- 2 



C« 2 3 



rBe 


Be 


a 


0 £ 


Be 


rBe 


0 


§ - 2 


r 


Be rBe 


T+ 1 

r 7 Ee 


- 

f^Be 



_o_ 

rBe 

a/9 

rBe 



2 a_ 
rBe 

-t. 

Be 

• 2/9 
r*Be 



M 



_ 2 § 

r 7 Ee 

T + 3 -a2 

r^He Be 



[M t ] 



1T-4ra+3fif_- 1 
r r JRe v 3 Re 

-aT- a 

if 7 Be 



ajB -4/9 
r*Be 



-aT- a - a/9 

r^Be r^Be 

_fl2 _ a z_ -5 t-_J3 

r^Be rBe r r^Be 

-£t+ 3/9_ !T + 2fi2_-_3 - 

r r^Ke r r^Be r^Be 



(D-55) 



(E-40) 



( E-4 1 ) 



(E-42) 



_a£ 

rBe 

(E-43) 
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[M o ] = 

^2 T -_4|2_ 
r 2 r*Ee 

J 

r r 2 Ee 



-aT+4ra2+ a 
r r^Ee 

tT+ a 2 -J3 2 
~le r*Ee 



r^Ee 



-afiT J? T + 4a£-_3fi-_2fit 
r r 2 r^fie r 2 Ee 



-a# T+ a/9 
r r^Ee 

-j3T-4a/3+ £ + 2a 2 £ 

r 2 r*Ee r 2 Ee 

(a 2 - 1 )T+ 3 + a 2 + 2/3* 

r 2 r T Ee r 2 Ee FEe 



(E-44) 



'V " 



t n, : 



CH 0 ) - 



1 0 


o“ 






0 0 


0 






0 0 


1 






1 ■ 
r 


-a 


o” 




-a 


0 


-£ 








r 




0 • 




1 






i 


r 




£| 


-a 


__ 




r 2 


r 




r 


0 


t 




A 








r 2 


“fifi 

r 


4 

r 2 


a 2 


-1 

r 2 



and 



T = 



aU - t 
Ee 



- 1 (« 2 ^£) 
Ee r 2 



2 a{1-r 2 ) 

and t is defined by eguation (E-34) . 



(E-45) 



( E— 46) 



(E-47) 



(E— 48) 
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APPENDIX F 



BOUNDARY CONDITIONS AT THE WALL 
1. Boundary Conditions for Plane Flow at y=±1 



The boundary condition is that the velocity v is 
zero at the walls. The boundary conditions for W must be 
determined. Ey the definition of W, 



V = VXW . 



(2-7) 



The form of v and W is given by equations (2-9a) r (2-9b) and 
(2-10). Equation (2-7) can be written 



l x \ 

u (y) e 




' Xn 

f (y) e 


, X 

( v (y) e 


\ = vx( 


x t 

g(y)e > 


X 

|w (y) e 




X 

h (y) e 

s / 




i 


j 


= 


3/3 x 


3/3 y 




X 

f (y) e 


g(y)e X 



k 

3/3 Z 

X 

h (y) e 



i 



3 k 



a 

f(y) 



D P 

g (y) h (y) 



x 

e 



= i(Dh-/3g)e + j(/3f-ah)e + k(ag-Df)e 



(F-1) 



Eguation (F-1) must 



be satisfied by each 



component 



(F-2a) 



separately. Dividing by e , the results are 
u (y) = Dh (y) - (y) 

v(y) = /3f(y) - ah (y) 

v(y) = ag(y) - Df(y) 

As described in Appendix C, one component of the 
vector potential W may be arbitrarily set egual to 

a. Case 1. . . f (y) =0 

At y = ±1, equations (F-2) become 
u (±1) = 0 = Dh (±1) - /3g (±1) 

v (±1) = 0 = -ah (±1) 

w (±1) = 0 = erg (±1) . 

These equations imply 
g(±1) = 0 

h (±1) = 0 

Dh (±1) = 0 . 

b. Case 2. . . g (y) =0 

At y = ±1, equations (F-2) become 
u(±1) = 0 = Dh (± 1) 

v (± 1 ) = 0 = /3f (±1) - ah(±1) 

w (±1) = 0 = -Df (±1) . 



(F-2b) 

(F-2c) 

velocity 

zero. 

(F-3a) 

(F-3b) 

(F-3c) 

(F-4a) 

(F-4b) 

(F-4c) 

(F-5a) 

(F-5b) 

(F-5c) 
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These equations imply 

/3f (± 1} = ah (±1) (F-6a) 

Df (±1) =0 (F-6b) 

Dh (± 1) = 0 . (F-6c) 

c. Case 3. . . h (y) =0 

At y = ±1, equations (F — 2 ) become 
u(±1) = -/%<±1) (F-7a) 

v(±1) = /3f(±1) (F-7b) 

w(±1) = ag(±1) -Df(±1) . (F-7c) 

These equations imply 

f (±1) =0 (F-8a) 

Df (±1) =0 (F-8b) 

g (±1) = 0 . (F-8c) 

2. Boundary Conditions for Pipe Flow at r = 1 

The boundary condition is that the velocity v is 
zero at the wall. The boundary conditions for K must be 
determined. By the definition of W, 

v = vxW . (2-7) 

The form of v and W is 
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v (x,r,0, t) 



(F-9) 



= (e u(r)+e v(r)+e w(r))e 
x r 0 



X 

W(x,r,€,t) = (e f(r)+e g(r)+e h (r) ) e 

x r 0 



where 



X = ax + /?0 + yt . 



In cylindrical coordinates, equation (2-7) can be 



/ Xx 


/ x \ 


u (r) e 




f (r) e 1 


X 




X 1 


<v (r) e 


> = «( 


9 (r) e 


X 


• 


X 


w(r)e 

\ i 




h (r) e 

' / 



1e 

r x 


1e 
r r 


e 

0 




1 e 

r x 


1 e 
r r 


e 

0 


3 

35 


HI 


d 


= 


a 


D 


P 


X 

f e 


X 

ge 


X 

rhe 




f 


g 


rh 



x x - 

= Jle (D(rh)-/?g)e + _le (^?f-arh) e + e^ 



Equation (F— 12) must be satisfied by each 

X 

separately. Dividing by e , the results are 



u (r) = JD (rh (r) ) - £g (r) 
r r 

= Jh(r) + Dh (r) - £g(r) 
r r 



v(r) = £f(r) -ah (r) 
r 



(F- 10) 

(F-11) 

written 



(ag-Df) 

(F- 12) 
component 

(F- 13a) 
(F-13b) 
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v (r) = ag (r) -Df (r) 



(F- 1 3c) 



As described in Appendix C, one component of 
arbitrarily set egual to zero. 

a. Case 1 . . . f (r) =0 

At r=1, eguations { F— 13) become 
u(1) = 0 = h(1) + Dh { 1 ) - (3g (1) 

v{1) = 0 = -ah { 1) 

w { 1) = 0 = ag (1) 

These eguations imply 
g{1) = 0 

h{1) = 0 

Dh { 1) = 0 



b. Case 2. . . g (r) =0 

At r=1/ eguations (F-13) become 
u (1) = 0 = h (1)' + Dh (1) 

v(1) = 0 = £f(1) -ah { 1) 

w ( 1) = 0 = -Df (1) 

These eguations imply 
Pf (1) = ah (1) 

Df ( 1) = 0 



may be 

(F- 14a) 
(F-14b) 
(F- 14c) 

(F-15a) 
(F- 15b) 
{F- 15c) 

(F- 16a) 
(F- 16b) 
(F- 16c) 

(F- 17a) 
(F- 17b) 
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h(1) = -Dh(1) 



(F- 17c) 



c. Case 3. . . h (r) =0 

At r=1, equations (F-13) become 
u (1) = 0 = -/% (1 ) 

v (1) = 0 = £f(1) 

w (1) = 0 = ag (1) -Df (1) 

These equations imply 
g(1) = 0 

f (1 ) = 0 

Df ( 1) = 0 



(F-18a) 
(F- 1 8b) 
(F- 1 8c) 

(F- 19a) 
(F- 1 9b) 
(F-19C) 
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APPENDIX G 



BOUNDARY CONDITIONS 
FOR PIPE FLOW AT r = 0. 

PERIODICITY OF ©. 

1. Boundary condition for 0 

The condition imposed is that functions be single 
valued at 0=0 and 0 = 27r. The assumed form of the solution is 



W(x,r,6,t) = [e f(r)+e g (r) +e h (r) ]e 

x r 0 



(E-1) 



where 

X = ax + Pe + yt . ( E— 2) 

Thus the 6 dependence may be separated as the factor 



pe 

e = cos (-i/36) + isin (-i/30) 



(G-1) 



This terra must be single valued at 0=0 and Q=2v. Both the 
real and imaginary parts must match. This results in two 
eguations. 



cos (-i/327 r) = 1 



(G-2a) 



sin (-i/327 r) = 0 



(G-2b) 



These conditions are met only for arguments of sine and 
cosine which are even multiples of rr, that is 



or 



(-i/327 r) = 2n7T, n=0,1,2,... 



$ - ni, n = 0, 1,2,... 



(G— 3) 



(G-4) 
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2. Single-Valued ness at r=0 



It is required that a vector function with the form 
of eguaticn (2-9) be single-valued at the origin, r=0. 

Consider a vector V which has the form 

n 



X 
e 

where 

X = ax + inG + y t . 



V = e u (r) + e v (r) + e w (r) 
n x n r n 6 n 



(G— 5) 



(G-6) 




Figure G-l Vector Function Rotated <p 
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Consider figure G-1 where point P is rotated an 
angle <J> about the origin to become point P 1 , both at a 
distance r from the origin. The unit coordinate vectors are 



<D 1 

N 

<t> 1 


and e at P . 


and 


e* / 


e • , and 


e * at 


P 1 . Vectors e 


x r 


0 




X 


r 


0 


X 


and e • 


extend out 


of 


the 


picture 


normal 


to e and e . At 


X 












r 0 


point P' 


, V becomes 














n 














1- 






_ 


\ X* 




V* 


= e' u (r) + 


e 1 


v (r) 


+ e* w 


(r) e 


(G-1) 


n 


l x n 


r 


n 


0 n 


1 





where 

X' = orx + in (0+$) + yt . (G— 8) 

Then 



X' X ini 

e = e e 



and 



(G-9) 



e * = e (G-lOa) 

x x 



e' = e cos<J) + e sin<}> (G-lOb) 

r r 0 



e 1 = e cos<J> - e sin<J> (G-IOc) 

0 0 r 



Substituting (G-9) and (G-10) into (G-7) and rearranging 
terms results in eguation { G— 11)- 



V' 

n 



e u (r) 
x n 



+ e (v (r) cos<J> 
r v n 



w 

n 



(r) sin<J>) 



+ e (v (r)sirnj) + 
© v n 



w (r) cos<t>) 
n ' 



X in<{i 
e e 



(G-11) 
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It is required that 



lira V = lira V' . (G-12) 

r-*0 n r-*0 n 



If equations (G- 1 1) and (G— 5 ) are set equal to each other at 
r=0, then each component of this vector equation must 
separately satisfy the equality. The result is the 
fcllovinq three equations. 



in<{> 

u (0) = u (0) e (G- 13) 

n n 



incj) 

v(0) = [v (0)cos<{> -w (0)sin<J>]e {G— 14) 

n n 



in<J> 

v (0) = [v (0)sin<l> + w (0)cos<J>]e {G— 15) 

n n n 

Equation (G-13) can be written 

in<J> 

u (0) [ 1-e ] = 0 . (G- 16) 

n 

in$ 

For n=0, 1-e is identically zero, so u q (0) is 

in<t> 

unrestricted. For n#0, 1-e is, in general, not equal to 

zero so u ( 0) = 0 for n*0. Equations (G-14) and (G-15) are 
n 

rewritten 



in<f> in<j> 

v (0) (1-cos<J>e ) + w (0) (sin<J>e ) = 0 
n n 



in<{> 

v (0) (-sin<Jie ) 
n 



in<}) 

+ w (0) (1-cos<$e ) 

n 



0 



(G- 17) 



(G- 18) 
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In (G— 1 7) and (G-18), v (0) and w (0) must equal zero unless 

n n 

the determinant of the coefficients vanishes. Setting the 
determinant of the coefficients equal to zero yields the 
following equation. 

in4> in4> 

(1-cos<J>e ) 2 + (sin4>e ) 2 = 0 (G-19) 

Multipling the terms out, using Eulers equation and noting 



that sin 2 <J> + ccs 2 $= 1, results in 

1 - 2cos$ (cosn<J>+isinn<>) + cos2n<t> + isin2n<}> = 0 (G— 20 ) 

The real and imaginary parts of equation (G-20) must 
separately te satisfied as follows. 

1 - 2cos(}>cosn4> + cos2n<J = 0 (G-21) 

-2cos<J>sinn$ + sin2n<}) = 0 (G-22) 



Using well known identities for cos2n<J> and sin2n<J> yields the 
following expressions. 

1 + cosn4> (cosn<j>-2cos4>) -sin 2 n<t> = 0 (G-23) 

sinn$ (cos<j>-cosn<}>) = 0 (G-24) 

For n=0, equation (G-2 h) is satisfied identically, but 



eguaticn (G-23) reduces to 

2 (1-cos$) = 0 (G-25) 

which, in general, is not satisfied for arbitrary <t>. For 

n=1, equations (G-23) and { G— 24 ) become 

1 - (sin 2 <$>+cos 2 4>) = 0 (G-26) 

sin$ (cos4)-cos4)) = 0 (G-27) 

which are satisfied identically for any value of <t>. For 



n>1, the equations will not be satisfied. In particular, 
equation (G — 2 4 ) requires that either sinn<{>=0 or cos4>=cosncJ> 
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which cannot be satisfied for arbitrary and for n>1. 



A relationship between v^ (0) and w^fO) is i 

the previous result. Using equation (G-17) to 
ratio of (0) and w^ (0) results in 

i<t> 

(0) sm^e sm(J> 

i<J> -i<J> 

w^ (0) l-cos^e e -cos(J» 



-sin<j> 

= i 

cos (-$) + isin (-<|>) -cos<J> 



or 



v (0) = iw (0) 

i i 

Equation (G — 1 8) yields the same result. A summary 
results of applying the restriction that a vector 
of the form specified by equations (G-5) and (G-6) 
valued at r=0 fellows. 

n=0 u (0) unrestricted 
o 

v (0) = 0 

o 

w (0) =0 

o 



n=1 u (0) = 0 

i 



v (0) = iw (0) 

i i 



mplied by 
find the 

(G- 28) 

(G-29) 

of the 
function 
be single 

(G-30a) 

(G-30b) 

(G-30c) 

(G-31a) 
(G — 31b) 
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(G-32a) 



n>1 u (0) = 0 

n 



v (0) = 0 

n 



(G-32b) 



w (0) = 0 

n 



(G-32c) 



3. Symmetry of Solution 



Consider figure G-1 with <t> = ir, 

(G-5) and (G— 11) 



From equations 



V = [e u (r) + e v (r) + e w (r) ]e 

n x n r n 0 n 



(G — 3 3) 



— - - - X in 7r 

V* = [e u (r) - e v (r) - e w (r) ]e e 
n x n r n 0 n 



(G-34) 



inir n - 

Note that e = (-1) , and let the components of V and V* 

n n 

be 0 , V , W and U', V' and W 1 . From equations (G-33) and 
n n n n n n 

(G-34) the components of V and V* have the following 

n n 

relations . 



U = 



O' (-1) 

n 



n 



(G-35) 



V = 
n 



-V' (-1) 
n 



n 



(G-36) 



n n 



(G-37) 



Thus the components of a vector with the form of V must act 

n 
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either as an even or odd function at the origin, depending 

upon n. But the functional dependence of U on r is 

n 

contained in u (r) . Thus u (r) must look like an even or 

n n 

odd function. Similarly, v and w must behave like even or 

n n 

odd functions. In summary, eguations (G-35) , (G-36) and 

(G-37) imply the following 



n even 


u (r) 
n 


even symmetry 


(G-38a) 




v (r) 1 
n 


> odd symmetry 


(G-38b) 




w (r) 
n ' 




(G-38c) 


n odd 


u (r) 
n 


odd symmetry 


(G-39a) 




v (r) 1 
n 


> even symmetry 


(G-39b) 




w (r) 




(G-39c) 



n 



4 . Summary 

If the curl operation is performed upon a vector 
with the conditions of symmetry from part 3, the resulting 
vector obeys the same rules of symmetry, that is, the 
specification that even or odd derivatives of the components 
equal zero is not changed by applying the curl. All vector 
functions are assumed to be of the form of equation ( G— 7) . 

The condition of single-valuedness at the origin for 
n=0 and n=1 has the same property. Using the boundary 
conditions at the origin from eguations (G-30) or (G-31) and 
applying the curl operator does not result in different 
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conditions on the results. This can be verified by 
performing the curl operation on a general vector of the 
form of (G-7) , imposing the boundary conditions on the 
original vector's components and observing the resulting 
restrictions on the components of the curl. 

For n>1, this is no longer true. Only the cases of 
n=0 and n=1 are considered in this paper. Thus, for n=0 and 
n=1, the vorticity, velocity and velocity vector potential 
will all have the boundary conditions previously discussed. 
The boundary conditions on the velocity vector potential at 
r=0 are summarized below. 

n=0 all odd derivatives of f = 0 
all even derivatives of g = 0 
(including g (0) = 0) 
all even derivatives of h = 0 

(including h (0) = 0) (G-40) 

n=1 all even derivatives of f = 0 
(including f (0) = 0) 
all odd derivatives of g = 0 
all odd derivatives of h = 0 

g (0) = ih (0) (G-4 1) 
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PROGRAM #1 



PROGRAM TO PRINT EIGENVALUES 
FOR THE 3-D POISEUILLE FLOW PROBLEM 



THIS PROGRAM SOLVES THE LINEARIZED NAVIER STOKES 
EQUATION FOR POISEUILLE FLOW. THE EIGENVALUES 
RESULTING FROM THE FINITE DIFFERENCE APPROXIMATION 
ARE PRINTEC. 

INPUT 

THE FOLLOWING MAY BE INPUT TO THE PROGRAM AS DATA 
USING NAMELIST, 'LIST'. NOTE, THE DEFAULT VALUES 
ARE ONLY SET INITIALLY AND VALUES SET BY THE 
USER WILL NOT BE CHANGED BETWEEN RUNS 

N - HALF OF THE NUMBER OF FINITE DIFFERENCE GRID 
POINTS ACROSS THE CHANNEL NOT INCLUDING THE END 
POINTS. N MUST BE .LE. MDIM, WHICH IS THE 
DIMENSION OF THE MATRICES IN THIS PROGRAM. 
DEFAULTED TO THE VALUE OF NDIM, THAT IS THE 
DIMENSION OF THE MATRICES. SEE PROGRAM BELOW FOR 
THE DEFAULT VALUE. 

REY - THE REYNOLDS NUMBER (REAL*8) DEFAULT 
VALUE = 60C0.0 

AR,AI,BR,BI - THE REAL AND IMAGINARY PARTS OF THE 
PERTURBATION WAVE NUMBERS (REAL*8) 

DEFAULTED TO O.O, 1.0, 0.0, AND 0.0 RESPECTIVELY 

VEL - THE VELOCITY OF THE MOVING COORDINATE 
REFERENCE SYSTEM FOR WHICH THE STABILITY IS 
DETERMINED. (REAL*3) DEFAULTED TO 0.0 

OUTPUT 

THE CUTPUT OF THIS PROGRAM IS A TABULATION OF THE 
EIGENVALUES. TWO LISTS ARE PRINTED, ONE FOR THE 
EIGENVALUES CORRESPONDING TO EVEN EIGENFUNCTIONS 
AND CNE FOR THOSE CORRESPONDING TO ODD EIGEN- 
FUNCTIONS. THE STABILITY OF EACH EIGENVALUE IS 
PRINTED AND THE LEAST STABLE E IGENVALUE IS MARKED 
WITH ASTERISKS. A PLOT OF THE EIGENVALUES IS ALSO 
PFINTED. 

SUBROUTINES 

THIS PROGRAM CALLS THE SUBROUTINE ' DEI GEO' TO 
SOLVE FOR THE EIGENVALUES. SUBROUTINE ’PLOTP* 

IS LSEO TO PLOT THE EIGENVALUES ON THE PRINTER. 



IMPLICIT REAL*8 <A-H,0-Z) 

CCMPLE X* 1 6 A , B 
REAL*4 GR4 ( C>0 ) , G I 4( bO ) 

REAL*8 GRE(60) ,GIE(60) ,GRO(60) , G I 0 ( 60 ) 
COMPLEX* 16 XMAT (60,30,3) 

C0MPLEX*16 YM AT (60,30), WVEC ( 60) ,BMAT(5,60) 
EQUIVALENCE ( Y MAT ( 1,1) ,XMAT( 1,1,3)), 

* (BMATI 1, 1) ,XMAT{1,1,3) ) , 

* ( WVEC( 1) , XMA T ( 1,6,3) ) 

NAMELIST / LIST / N , RE Y , AR , AI , BR , B I , VE L 

INITIALIZE VARIABLES (SET DEFAULT VALUES) 
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MDIM = 60 
N = 60 
REY = 6CG0D0 
AR = 0C0 
A I = ICC 

ep = oco 

B I = 0D0 

READ NAMELIST AND SET ALPHA AND BETA 

1 READ(5, LIST, END-100) 

A = DCMFLXtAR, AI ) 

B = DCMPLX(BR,BI ) 

PRINT INPUT VALUES AS PAGE HEADING FOR EIGENVALUE LIST 

WRITE(6,9004) 

9004 FORMAT (' 1' ) 

WRI T E ( 6 » 9005 ) N , PE Y , A , B , VE L 

9005 FORMAT ( * N = ',I4,/,* REY = • ,F10 . 2, 8X, • ALPHA =• , 

* 2F12. 7,8X, 'BETA = • , 2F 1 2. 7 , / , ' VEL =',F7.2) 

CALL SUBROUTINE TO SOLVE FOR EIGENVALUES. 

CALL DEIGEOI A, B, REY, N, MDIM ,GRE ,GIE ,GRO,G 10 ,XMAT , YMAT, 

* BMAT , WVEC) 

DETERMINE WHICH EIGENVALUE 

TEMP = - 1D10 
MARK = 1 



IS THE LEAST STABLE. 



DO 20 1=1, N 

IF (GRCU )+AR*VEL.LT.TEMP ) GO TO 
TEMP = GRO(I )+AR*VEL 
I TEMP = I 
20 CONTINUE 



20 



DO 40 1=1, N 

IF ( GRE ( I ) +AR^VEL .LT.TEMP) GO TO 40 
TEMP = GRE (I )+AR*VEL 
ITEMP = I 
MARK = 2 
40 CONTINUE 

LIST EIGENVALUES FOR ODD EIGENFUNCTIONS 
WRITE { 6,9003) 

9003 F0RMAT(///,6X, 'GAMMA REAL* , 5X, ' GAMMA I MAG* , 1 2X , ' ST A3 • ) 
WR I T E ( 6 , 9C06 ) 

9006 FORMAT ( ' OE IGE N VALUES FOR ODD EIGENVECTORS',/) 

CO 50 1=1, N 

TEMP = GPO( I )+AR*VEL 

WRITE (6, 9000) GRO ( I ) , G I 0 ( I ) , TEMP 

I F ( I .EC. I TEMP. AND. MARK. EQ. 1) WR I TE ( 6 , 900 1 ) 

5C CONTINUE 

9000 FORMAT ( • 0 ' , 1P2D1 5. 4, 1PD20.4) 

9001 FORMAT (•+• ,52X, •***• ) 

LIST EIGENVALUES FOR EVEN EIGENVECTORS. 

WRITE! 6,9007) 

9007 FORMAT ( • 0 E IGENVALU ES FOR EVEN EIGENVECTORS',/) 

DO 55 1=1, N 

TEMP = GPE< I )+AR*VEL 

WR I T E ( 6 , 9000 ) GRE ( I ) , G I E ( I ) , TEMP 

IF ( I . EC. ITEMP. AND. MARK. EQ. 2) WRITE (6 ,9 0" l ) 

55 CONTINUE 

PUT EIGENVALUES INTO SINGLE PRECISION VECTORS TO PASS TO 
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SUeRCUTINE TO DO PLOTTING FOR ODD FUNCTIONS. 

DO 60 1 = 1. N 

GR4 ( I ) = SNGL ( GRO ( I ) ) 

6C GI4( I) = SNG L ( G I 0( II ) 

WRITE (6,9004) 

CALL FLCTP(GR4»GI4,N»0) 

WRITE( 6,9005) N , RE Y , A , B , VE L 

SIMILARLY PLOT EIGENVALUES FOR EVEN EIGENFUNCTIONS. 

DO 6 5 I = 1 , N 

GR4 ( I) = SNGL (GRE( I ) ) 

65 G 14 ( I) = SNGL ( G I E ( I ) ) 

WRITE( 6,9004) 

CALL PLCTP ( GR4 »G 14 , N , 0 ) 

WRITE( 6,9005) N, RE Y, A, B, VEL 

GO TO 1 

IOC WRITE( 6,9004) 

STOP 

END 
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PRCGRAM #1A 

PROGRAM TO PRINT EIGENVALUES 
FOR THE 3-D POISEUILLE FLOW PROBLEM 
AND PLOT CN CALCCMP PLOTTER 



THIS PFCGRAM SOLVES THE LINEARIZED 
EGUATICN FOR POI SEUSLLE FLOW. THE 
RESULTING FROM THE FINITE DIFFEREN 
ARE PRINTED. 

INPUT 

THE FOLLOWING MUST BE INPUT TO 
USING NAMELIST, 'LIST' . 

N - HALF OF THE NUMBER OF FINI 
POINTS ACROSS THE CHANNEL NOT 
POINTS. N MUST BE ,LE. MDIM, 
DIMENSION OF THE MATRICES IN T 
DEFAULTED TO THE VALUE OF MDIM 
DIMENSION OF THE MATRICES. SE 
THE DEFAULT VALUES. 

REY - THE REYNOLDS NUMBER (RE 
DEFAULT VALUE = 6000. 

A R » A I , BR , B I - THE REAL AND IMA 
PERTURBATION WAVE NUMBERS (PE 
TO 0., 1., 0. AND 0. REoPECTIV 

VEL - THE VELOCITY OF THE MCVI 
REFERENCE SYSTEM FOR WHICH THE 
DETERMINED. (REAL*8). DEFAUL 

GR SCL , GI SC L - THE X AND Y SCAL 
GRAPH IN UNITS PER INCH. BOTH 
SEE SUBROUTINE DRAW FOR REST RI 
VALUES 

NW ,NH - THE WIDTH AND HEIGHT 0 
IN INCHES. DEFAULTED TO 6 AND 
SEE SUBROUTINE DRAW FOR RESTRI 
VALUES 

TITLE - A VECTOR OF LENGTH 12 
WITH THE GRAPH (REAL*8) 

OUTPUT 

THE OUTPUT OF THIS PROGRAM IS 
EIGENVALUES. TWO LISTS ARE PR 
EIGENVALUES CORRESPONDING TO E 
AND CNE FOR THOSE CORRESPONDIN 
FUNCTIONS. THE STABILITY OF E 
PRINTED AND THE LEAST STABLE E 
WITH ASTERISKS. THE E IGENVALU 
THE CALCCMP PLOTTER. 

SUBROUTINES 

THIS PROGRAM CALLS THE SUBROUT 
SOLVE FOR THE EIGENVALUES AND 
TO GRAPH THE RESULTS. 



NAVIER STOKES 
EIGENVALUES 
CE APPROXIMATION 



THE PROGRAM AS DATA 



TE DIFFERENCE GRID 
INCLUDING THE END 
WHICH IS THE 
HIS PROGRAM. 

, THAT IS, THE 
E PROGRAM BELOW FOR 



AL*8 ) 



GINARY PARTS OF THE 
AL*8). DEFAULTED 
ELY. 

NG COORDINATE 
STABILITY IS 
TED TO 0. 

E OF THE OUTPUT 
DEFAULTED TO 0.2. 
CTIONS ON THEIR 



F THE OUTPUT GRAPH 
8 RESPECTIVELY. 
CTIONS ON THEIR 



WHICH IS PRINTED 



A TABULATION OF THE 
INTED, ONE FOR THE 
VEN EIGENFUNCTIONS 
G TO ODD EIGEN- 
ACH EIGENVALUE IS 
IGENVALUE IS MARKED 
ES ARE GRAPHED BY 



INE ' D El GEO' TO 
SUBROUTINE 'DRAW* 
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IMPLICIT REAL*8 (A-H,0-Z) 

CGMPLEX*16 A , B 
REAL*4 GR4(30) ,GI4(30) 

REAL* 8 GRE(60),GIE(60) ,GR0(60) ,GI0(60) 

COMPLEX* 16 XMAT( 60,30, 3) 

CC MPL EX*1 6 YMAT (60,30) , WV EC ( 60) ,8MAT( 5,60) 
EQUIVALENCE (YMAT( 1 , I) ,XMAT(1, 1,3) ), 

* (BMAT( 1,1) ,XMAT( 1,1,3) ) , 

* ( WVEC( I) ,XMAT( 1,6,3) ) 

RE AL*8 LABEL, TITLE ( 12 J 

REAL*4 GR SCL , G I SCL 
DATA LABEL/* •/ 

CATA TITLE/' HARRISON* , ' 0279 ',10*' •/ 

NAMELIST / LIST / N , REY , AR , A I , BR ,8 I , VEL , GRSCL , GI SCL , 

* NWjNH, TITLE 

INITIALIZE VARIABLES (SET DEFAULT VALUES) 

MDIM = 60 
N = 60 

REY = 6C00D0 
AR = OCO 
AI = ICO 
BR = 0D0 
B I = OCO 
VEL = ODO 
GRSCL = .20 
G I SCL = .20 
NW = 6 
NH = 8 

READ NAMELIST AND SET ALPHA AND BETA 

1 READ(5,LIST,END=I00) 

A = DCF FLX ( AR , AI ) 

B = DCMPLX (BR, BI ) 

PRINT INPUT VALUES AS PAGE HEADING FOR EIGENVALUE LIST. 

WRITE ( 6 ,9004) 

9004 FORMAT Cl') 

WRI TE ( 6 , 9CC5 ) N , RE Y , A, B , VE L 

9005 FORMAT ( * N =*,I4,/,' REY = • , F 10 . 2, 8X , • ALPH A =* , 

* 2F12.7,8X, 'BETA = * , 2F1 2 . 7 , / , • VEL =*,F7.2) 

CALL SL3FCLTINE TO SOLVE FOR EIGENVALUES 

CALL DEIGEG( A, B, REY, N, MDIM , GRE , GI E ,GRO , G 10 ,XMAT , YMAT, 

* BMAT , WVEC ) 

DETERMINE LEAST STABLE EIGENVALUE 

TEMP = - 1 DIO 
MARK = 1 

CO 20 1=1, N 

IF (GFO( I )+AR*VEL.LT. TEMP) GO TO 20 
TEMP = GRO(I )*AR*VEL 
ITEMP = I 
20 CONTINUE 

CO 40 1=1, N 

IF (GRE ( I ) +AR*VEL .LT.TEMP) GO TO 40 
TEMP = GRE ( I ) +AR*VEL 
ITEMP = I 
MARK = 2 
40 CONTINUE 

LIST EIGENVALUES FOR ODD EIGENVECTORS 
WR I T E ( 6 , 9003 ) 

9003 FCRMAT(///,6X, 'GAMMA REAL 1 > 5X f 1 GAMMA I MAG* , 1 2X , ' ST AB • ) 
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WR ITE ! 6 , 9006 ) 

90C6 FORMAT ( 'OEIGEN VALUES FOR ODD EIGENVECTORS',/) 

C 

CO 50 1=1, N 

TEMP = G PO ( I )+AR*VEl 

WRITE! 6,9000) GRC ( I ) , G I 0 ( I ) , TEMP 

IF (I .EC. ITEM P. AND. MARK. EQ.l) WR I TE ! 6 , 900 1 ) 

50 CCNTINIE 

90C0 FORMAT! 'O' , 1P2D1 5. 4, 1PD20. 4) 

9001 FORMAT (•+• ,52X, '***' ) 

LIST EIGENVALUES FOR EVEN EIGENVECTORS. 

WRITE! 6,9007) 

90C7 FORMAT! • OE IGENVALUES FOR EVEN EIGENVECTORS',/) 

DC 55 1=1, N 

TEMP = GRE(I )+AR*VEL 

WRITE! 6,9000) GR E ( I ) , G I E ( I ) , T EMP 

IF (I . EC. I TEMP. AND. MARK. EQ. 2) WR I TE ( 6 , 900 1 ) 

55 CONTINUE 
C 

C PUT EIGENVALUES INTO REAL*4 VECTOR IN GROUPS OF 30 OR LESS 
C TO PASS TO SUBROUTINE DRAW. CALL DRAW REPEATEDLY UNTIL 
C ALL EIGENVALUES FOR ODD EIGENVECTORS ARE PLOTTED. DO THE 
C SAME FOR TFE EIGENVALUES FOR EVEN EIGENVECTORS. 

C 

NGCOCE = 1 
IE = 30 
I MODE = 0 

IF(30.GT.N) IE = N 
C 

60 DO 65 1=1, IE 

GR4 ! I ) = SNGLIGRC! IM0DE*30 + I)) 

65 G 1 4! I ) = SNGL(GIC( IM0DE-30 + I) ) 

C 

CALL CRAW! I E , GR4 , G 14 , NGC ODE , 3, LABEL, TITLE, 

* GRSCL , GI SCL , NH >NW , 2 , 2 , NW , NH , 0, LAST ) 

C 

NGCODE = 2 

IF t ( IMCDE*1)*30+1.GT.N) NGCODE = 3 
DO 7 0 1 = 1, IE 

GR4 ( I ) = SNGL! GRE! 1M0DE*30*I ) ) 

7C G I 4! I ) = SNGLtGIE! IMCDE*30 + I) ) 

C 

CALL CRAW! I E , GR4 , G 14 , N GCODE , 1 , L A BEL , T I TL E , 

* GRSCL,GISCL,NH,NW,2,2,NW,Nh,0, LAST) 

C 

IF ( NGCCCE. EQ.3 ) GO TO 80 
I MODE = I MODE + 1 

IF! ( IMCDE+1)*30.GT.N ) IE = N - IM0DE*30 
GO TO 60 
C 

80 GO TO 1 
100 STOP 
END 
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PROGRAM HZ 

PROGRAM TO FIND THE SHAPE OF 
THE STABILITY CURVES 
FOR 3-D POISEUILLE FLOW 



THIS PROGRAM FINDS THE SHAPE OF THE STABILITY CURVES. 
IT PERFORMS FIVE FUNCTIONS DEPENDING ON THE VALUE OF 
THE INPUT PARAMETER ‘MODE'. INPUT IS BY THE 
NAMELIST ’LIST*. OTHER INPUT AND OUTPUT ARE 
DETERMINED BY THE VALUE OF MODE AS WELL AS THE 
FUNCTION PERFORMED. IN ALL CASES, ONLY THE PARAMETERS 
A I AND REY ARE VARIED. OTHER PARAMETERS ARE ALWAYS 
hELD CONSTANT THROUGHOUT THE CALCULATIONS. 

MODE = 1 

THE STABILITY AT ANY GIVEN POINT IS FOUND 

INPUT 

AR,AI,BR,BI - THE REAL AND IMAGINARY PARTS 
OF ALPHA AND BETA. (REAL*8) 

VEL - THE VELOCITY OF THE MOVING COORDINATE 
FRAME FOR WHICH THE STABILITY IS TO BE 
DETERMINED. <REAL*8) 

REY - THE REYNOLDS NUMBER. (REAL*8) 

OUTPUT 

THE PROGRAM PRINTS THE INPUT AND THE RESULTING 
STABILITY. 

MODE = 2 

FOR A GIVEN VALUE OF REY, AI IS ADJUSTED TO FIND THE 
POINT WHICH HAS THE DESIRED STABILITY. 

INPUT 

AR,BR,BI,REY ,VEL - VALUES TO BE HELD CONSTANT 
DURING THE SEARCH FOR THE VALUE OF AI. 

AI - AN INITIAL GUESS. 

CAI - AN ESTIMATE OF THE ERROR IN THE GUESS 
FOR AI. ( R EAL* 8 ) 

STAB - THE STABILITY FOR WHICH THE VALUE OF AI IS 
DESIRED. (REAL-3) 

OUTPUT 

THE INPUT VALUES INCLUDING THE INITIAL GUESS FOR 
AI ARE PRINTED ALONG WITH THE FINAL APPROXIMATION 
FOR AI. 

MODE = 3 

FOR A GIVEN VALUE OF AI, REY IS ADJUSTED TO FIND THE 
POINT WHICH HAS THE DESIRED STABILITY. 

INPUT 

AR , A I , BR , B I - VALUES TO BE HELD CONSTANT DURING 
ThE SEARCH FOR THE VALUE OF REY. 
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REY - AN INITIAL GUESS. 

DREY - AN ESTIMATE OF THE ERROR IN THE GUESS FOR 
REY. IREAL*8) 

STAB - THE STABILITY FOR WHICH THE VALUE OF REY 
IS CESIRED. 

OUTPUT 

THE INPUT VALUES INCLUDING THE INITIAL GUESS 
FOR REY ARE PRINTED ALONG WITH THE FINAL 
APPROXIMATION FOR REY. 

MODE = 4 

TFE POINT CF MINIMUM REY ON THE NEUTRAL STABILITY 
CURVE IS DETERMINED. 

INPUT 

AP,BR,BI,VEL - VALUES TO BE HELD CONSTANT DURING 
THE DETERMINATION OF THE MINIMUM REY. 

A I » REY - INITIAL GUESS FOR THE POINT OF 
MINIMUM REY. 

CA I , D REY - AN ESTIMATE OF THE ERROR IN THE 
GUESSES FOR AI AND REY. 

OUTPUT 

ThE PROGRAM PRINTS THE INITIAL GUESS, THE 
SUCCESSIVE APPROXIMATIONS AND THE FINAL 
ESTIMATE FOR THE VALUES OF AI AND REY. 

MODE = 5 

THE POINT CF MAXIMUM INSTABILITY IS DETERMINED. 

INPUT 

AR i BR , BI , VEL - VALUES TO BE HELD CONSTANT 

AI , REY - INITIAL GUESS FOR THE POINT OF MAXIMUM 
INSTABILITY. 

CAI , DREY - AN ESTIMATE OF THE ERROR IN THE 
GUESSES FOR AI AND REY. 

OUTPUT 

THE PROGRAM PRINTS THE INITIAL GUESS, THE 
SUCCESSIVE APPROXIMATIONS AND THE FINAL ESTIMATE 
FOR ThE VALUES OF AI AND REY. 



SUBROUTINES REQUIRED 

THIS PROGRAM CALLS THE SUBROUTINE ' D FI T2 ' AND USES 
THE FUNCTION ' El GLS* . 



IMPLICIT REAL*8 (A-H,0-ZJ 

NAMELIST / LIST / RE Y , AR , A I , BI , DAI ,D RE Y , VEL, BR, MODE , 
* STAB 
WFITE(6,9000) 

9000 FORMAT ( • 1* ) 

C 
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INITIALIZE VARIABLES 

STAB = CDO 
VEL=0D0 
BR = 0C0 
B I = 0C0 
AR = 0C0 
A I = ICO 
C A I = .010C0 
DREY = 250D0 
MCDE = 1 

READ NAMELIST VARIABLES 

1 READ(5 ,LIST, END=100) 

PEY1 = REY 
AITEMP = A I 

JUMP TC TYPE OF SOLUTION CALLED FOR BY VALUE OF MODE 
GC TO 110,20,30,40 ,50) , MODE 

MODE = 1 

FIND STAB FOR SINGLE POINT 

10 WRITE (6,9100) REY , AR, AI , BR , BI , VEL 
9100 FORMAT!///, 21X,' STABILITY AT POI NT ...',/, 2 IX ,' REY = ', 

* F8 . 1 , / ,2 IX , * AR =• ,F12.5,5X,'AI = ' , F 1 2. 5, / , 21X , 

* 'BR =' »F 1 2 . 5 » 5X , 'BI = ' ,F 12. 5 ,/ , 2 IX , ' FOR VEL =', 

* F6.2) 

USE EIGLS TO FIND STABILITY 

G 1 = EIGLS(REY,AR, AI ,BR,BI ,VEL) 

WR I TE l 6 , 9 101 ) G 1 

101 FORMAT !/,21X,' STAB =',1PD12.4) 

GC TO 1 



MODE = 2 

FIND AI GIVEN REY AND STAB 



2 C WR IT E ( 6 ,9200 ) STAB , REY , A I , BR , B I , AR , VEL , D AI 
920C FORMAT (///,21X»'FOR STAB = • , F9 . 5 , 5 X , • AT REY =',F8.1, 

* /,21X, 'INITIAL GUESS.. . AI = ' ,F 12 .5 ,/ ,21X , 

* 'GIVEN BR =' »F 12 .5 , 5X, 'BI = ' , F 12 . 5 , / , 2 7X , ' AR =', 

* F12.5,5X, 'VEL =' ,F6.2, / , 21X, 'DAI =',F10.5,/) 

CALCULATE STABILITY FOR INITIAL GUESS. 

G1 = EIGLS(REY,AR, AI ,BR,BI ,VEL) 

WR ITE ( 6 , 5201 ) A I , G 1 

9201 FORMAT ( 2 IX ,' FOR AI = • , F9 .5 , 5X , ' STAB =',1PD12.4) 

ADD C A I TO INITIAL AI AND DETERMINE STABILITY. 

All = AI + DAI 

G2 = EIGLS(REY ,AR, AI 1 , BR,BI ,VEL) 

WRITE(6»S201 ) A 1 1 , G2 

USE LINEAR APPROXIMATION TO FIND NEW VALUE OF AI . 
CALCULATE STABILITY. 

AI 2 = All + ( AI1-AI ) *( STAB-G2)/(G2-G1) 

G3 = EIGLS(REY ,AR,AI 2,BR,BI ,VEL) 

WRITE(6,5201) AI2,G3 

USE SUBROUTINE DFIT2 TO OBTAIN SECOND ORDER APPROXIMATION 
FOR FINAL VALUE OF AI. 

CALL DFIT2(AI ,AI 1, AI2,G1 ,G2,G3, STAB, AIF, DUMMY , DUMY 2 ) 
WRITE! 6,92C2) A I F 

9202 FORMAT (/ ,21X ,' FINAL ESTIMATE... AI =',F12.5) 
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C MODE = 3 

C FIND REY GIVEN AI AND STAB 

C 

30 WRITE(6,9300) STAB , A I , RE Y , BR , B I , AR , VEL , DREY 
9300 FORMAT (/// ,21X, • FOR STAB = • , F9. 5 ,5X, • AT AI =',F12.5, 

* /»21X> a INITIAL GUESS... REY = • , F8 . 1 , / , 2 IX , 

* 'GIVEN BR = ' ,F12.5,5X, '31 = ' ,F 12 . 5 , / ,27X , 

* 'AR = * , F 1 2 . 5 , 5X , 'VEL = ' , F6 .2 ,/ ,2 IX , ' DREY = ', 

* F8.lt/) 

C 



C CALCULATE STABILITY FOR INITIAL GUESS 
C 

G1 = E IGLStREY ,AR,AI ,BR, BI ,VEL) 

WRITE (6,9301 ) REY,G1 

9301 FGFMAT (21X,'F0R REY = • , F 8 . 1 , 5X , ' ST AB =',1PD12.4) 

C 

C ACC CREY TC INITIAL REY AND DETERMINE STABILITY. 

C 

REY 1 = REY + DREY 

G2 = EIGLS(REY1,AR,AI,BR,BI,VEL) 

WRITE (6,9301 ) RE Y1 , G2 
C 

C USE LINEAR APPROXIMATION TO FIND NEW VALUE OF REY. 

C CALCULATE STABILITY. 

C 

REY2 = REY1 + ( R EY 1 -R E Y ) * ( STAB-G2 ) / ( G2~G 1 ) 

G3 = EIGLS(REY2,AR,AI,BR,BI,VEL) 

WRITE(6,9301 ) REY2 ,G3 
C 

C USE SUBROUTINE DFIT2 TO OBTAIN SECOND ORDER APPROXIMATION 
C FOR FINAL VALUE OF REY. 

C 

CALL DFIT2(REY,REY1,REY2 ,G1 ,G2,G3, ST AB , REYF , DUMMY , 

. * CUMY2) 

WRIT E ( 6 ,9302 ) REYF 

9302 F0RMAT(/,21X, 'FINAL ESTIMATE... REY =',F8.1) 

GO TO 1 

C 



C MODE = 4 

C FINC POINT OF MINIMUM REYNOLDS NUMBER 

C ON THE NEUTRAL STABILITY CURVE 

C 



C OUTPUT VALLES FOR INITIAL GUESS 
C 

40 WR I T E ( 6 , 9400 ) A I , R EY , D A I ,D RE Y , B I , A R , VE L , BR 

9400 FCRMAT( • 1* ,////, 21X, 'POINT UF MINIMUM REY ON NEUTRAL', 

* ' STABILITY CURVE' ,/ ,2 IX, 

* 'INITIAL GUESS... AI =',F12.5,5X, 

* ' R EY 1 =' ,F8.1,/»21X,'DAI = • , FI 2. 5 , 5X , • DREY =• , 

* F8.1 ,/,21X, • BI =' ,F12.5,5X, ' AR =',F12.5, 

* / ,21X, 'VEL =' ,F6.2,5X, 'BR =' ,F 12.5 ,/// } 

C 

C FOR INITIAL REYNOLDS NUMBER, CHOOSE TW C MORE VALUES FOR 
C AI AND USE PARABOLIC CURVE FIT TO FIND MINIMUM STABILITY 
C WITH CORRESPONDING VALUE FOR AI 
C 

AIM = AI-CAI 
AI P = AI+DAI 

G1 = E IGLS(REY1, AR ,AIM,BR, BI ,VEL) 

WRITE ( 6,9401 ) AIM,G1 

9401 FORMAT (21X,*REY1.. • AI = ' , F 12. 5 , 5X , ' FOR STAB =', 

* 1PD12.4) 

G2 = EIGLS( REY1 ,AR,AI, BR,BI ,VEL) 

WRITE( 6,9401 ) AI,G2 

G3 = EIGLS(REY1,AR,AIP,BR,BI ,VEL) 

WRIT E ( 6 ,9401 ) AIP,G3 

CALL DFIT2( AIM, A I , A IP ,G1 ,G2,G3 , 1D2 , DUMMY ,AI 1 ,GMIN1 ) 

WRI TE ( 6 , 9402 > AI1.GMIN1 

9402 FORMAT (• 0 • ,20X, 'AI = • , F 1 2. 5 , 5X , ' FOR MAX G =• , 

* 1PC15.4, /) 
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CHCCSE NEXT GUESS FOR REYNOLD'S NUMBER USING DREY 
REY2 = REY1+DREY 

I F (GMI N1 .GT.ODO) REY2 = REY1-DREY 
AI = All 

WR ITE ( 6 , 9403 ) A I » REY2 

9403 FORMAT (///,21X,‘ SECOND GUESS... AI =',F12.5,5X, 

* ' R EY 2 =' » F 8 . 1 » / ) 

FIND MINIMUM STABILITY AND CORRESPONDING AI FOR SECOND 
GUESS OF REY 

AIM = AI-DAI 
AI F = AI+CAI 

G1 = EIGLSIREY2, AR , A I M , BR , B I , VEL ) 

WRITE (6, 9404) AIM,G1 

9404 FORMAT ( 21X » ’ REY2.. . AI = ' , FI 2. 5 , 5X , • FOR STAB = ', 

* 1PD12.4) 

G2 = EIGLS(REY2,AR, AI,BR,BI ,VEL) 

WRITE( 6,9404) AI,G2 

G3 = EIGLS(REY2fAR,AIP,BR,BI,VEL) 

WRITE(6,S404) AIP,G3 

CALL DFIT2(AIM,AI ,AIP,G1,G2,G3 ,1D2, DUMMY ,AI2 ,GMIN2) 

WRI TE(6,9402) AI2,GMIN2 

USE LINEAR FIT TO DETERMINE THIRD GUESS FOR REY AND AI 

REY3 = REYl-{ REY1-REY2 ) *GMI Nl/ ( GMI Nl-GMI N2 ) 

AI = AI1-(AU-AI2)*GMIN1/<GMIN1-GMIN2) 

WRITE (6,9405 ) A I , R EY3 

9405 FORMAT (///,21X, ' THIRD GUESS... AI =',F12.5,5X, 

* 'REY3 =• tF8.1>/) 

REPEAT PROCEDURE TO FIND MINIMUM STABILITY AT THIRD GUESS 
FOP PEY 

AIM = AI-DAI 
AI P = AI+CAI 

G1 = EIGLS(REY3,AR,AIM,BR, BI»VEL) 

WPITE(6,9406) AIM,G1 

9406 FORMAT ( 2 IX t * REY3 . . . AI = • , F 12. 5, 5X , ' FOR STAB =' , 

* 1PC12.4) 

G2 = EIGLS(REY3, AR , AI , BR,BI ,VEL) 

WRITEt 6,9406) AI,G2 

G3 = EIGLS(REY3,AR,AIP,BR,BI ,VEL) 

W RIT E ( 6 ,9 406 ) AIP,G3 

CALL DFIT2(AIM,AI ,AIP,Gl,G2, G3 , 1 D2 , DUMMY , A 1 3 , GMI N3 ) 

WRI TE ( 6 ,9402 ) AI3,GMIN3 

USE PARAEGLIC FIT TO DETERMINE FINAL GUESS FOR REY AND AI 

CALL CFIT2(AI1 ,A 1 2, A 1 3 , G MI N 1 , GM I N2 ,GMIN3 ,0D0,AIF, 

=* DUMMY , CUMY2 ) 

CALL DFI T2 ( REY 1 , REY2 , REY3 , GMIM 1 » GMIN2, GM IN3, ODO, REYF, 

* DUMMY , DUMY 2 ) 

WRIT E ( 6 , 9407 ) A I F, REYF 

94C 7 FORMAT ( /// , 2 IX , * FINAL ESTIMATE... AI =',F12.5,5X, 

* 'REY = s ,F8.1 ,////////) 

AI = AITEMF 

GO TO 1 

MODE = 5 

FIND POINT OF MINIMUM STABILITY 
(THAT IS, FOR MAX VALUE OF STAB) 

50 WRITE ( 6,9500 ) A I , REY , D AI ,D REY , B I , A R, VE L , BR 
950C FORMAT (• 1 ',///, 20X ,' POINT OF MINIMUM STABI LI TY' , / ,2 IX , 

* 'INITIAL GUESS... AI = ' , F 1 2 .5 , 5X , • LEY 1 =', F3.lt 

* / , 2 1 X , 'DAI =' ,FL2.5,5X, 'DREY = ' , F8 . L ,/,21X, 

* 'B1 = ‘ ,F12.5,5X, 'AR = • , F 12 . 5 , / , 2 IX , * VEL =',F6.2, 

* 5X , 'BR =• ,FI2.5,///) 
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c 

AIM = AI - DAI 
AIP = AI + DAI 

G1 = EIGLS(P.EY1,AR »AIM»BR»BI » V E L ) 

WRITE! 6,S5C1 ) A I M , G 1 

9501 FCFMAT ( 2 1 X , ' R E Y 1 • . • AI = ' ,F1 2 .5 , 5X, • FOR STAB = ', 

* 1PD12.4) 

G 2 = EIGLS(REY1,AR,AI,BR,BI,VEL) 

WRITE( 6,9501 ) AI,G2 

G3 = EIGLS(REY1, AR,AIP ,BR,BI ,VEL) 

WRITE(6,95C1) A I F , G3 

CALL DFIT2(AIM,AI ,AIP,G1,G2,G3 » 1D2 , DUMMY , A 1 1 ,GMIN1) 
WRITE (6, 9502 ) AI1.GMIN1 

9502 FORMAT ( * O' »20X, * AI = ' , FI 2. 5 , 5X , ' FOP STAB = ' , 1 PD1 5 .4 , / ) 
C 

REY2 = REY1 + DREY 
AI = All 

WRITE(6,9503> AI.REY2 

9503 FCRMAT (/// ,21X , ' SECOND GUESS... AI = ',F12.5,5X, 

* 1 REY 2 =' ,F8. 1,/) 

C 

AIM = AI - DAI 
AIP = AI + DAI 

G1 = EIGLS(REY2,AR»AIM,BR,6I,VEL) 

WR I T E ( 6 » 9 504 ) AIM,G1 

9504 FCFMAT (21X, 1 REY2.. . AI = ' , F 1 2 .5 , 5X , • FOR STAB =', 

* 1 PD 12 . 41 

G2 = EIGLS(REY2»AR,AI,BR»BI , VE L ) 

WRITE ( 6 ,95041 AI,G2 

G3 = EIGLSIREY2, AR,AIP ,BR,BI ,VEL) 

WRI TE ( 6 , 9504 J AIF,G3 

CALL 0FIT2(AIM f AI , AIP,G1 ,G2,G3 , 1D2 , DUMMY ,AI2 ,GMIN2) 
WRITE (6, 9502) AI2,GMIN2 
C 

REY3 = REY 1 - DREY 

IF (GMIK2 .GT.GMIN1) REY 3 = REY2 + DREY 

AI = A 1 2 + (AI1-AI2)*( REY3-REY2)/( REY1-REY2) 

WR I TE ( 6 , 9505 1 A I , R EY 3 

95C5 FORMAT!///, 21X, 'THIRD GUESS... AI =',F12.5,5X, 

* * REY 3 =' , F 8 . 1 , / ) 

C 

AIM = AI - DAI 
AIR = AI + DAI 

Cl = EIGLS(REY3, AR ,AIM,BR, BI ,VEL) 

WR IT E ( 6 , 9506 ) AIM,G1 

9506 FORMAT ( 2 IX , 1 REY3.. . AI = • , FI 2 .5 , 5X, ' FOR STAB =', 

* 1PD12.4) 

G2 = E I GL S ( RE Y 3 , AR , A I » BR , 3 I , VEL ) 

WRI T E ( 6 , 9 506 ) AI,G2 

G3 = EIC-LS (REY3,AR ,AIP , BR , B I ,VEL) 

WRITE(6,9506) AIP,G3 

CALL DFIT2IAIM, A I ,AIP,G1 ,G2,G3 ,1D2 , DUMMY ,AI3 »GMI M3) 
WRIT5(6,S502) AI3,GMIN3 
C 

CALL DFIT2(REY1,REY2,REY3,GMIN1,GMIN2,GMIN3, IDO, DUMMY , 

* REYF , GF ) 

C 

A = ( (All — AI2J/(REY1 -REY 2) -( AI 1-AI 3) / ( REY1-REY3) ) 

* /( (REY1**2-REY 2**2) /(REY 1-RE Y2 ) 

* - ( REY 1**2- RE Y3** 2) /( REY1-REY3) ) 

B = ( A 1 1 - A 1 2 - A* ( RE Y 1 ** 2-RE Y2**2 ) )/(REYl-REY2) 

C = AI 1- A*REY1**2-B*REY1 
AI F = A*REYF**2+B*REYF+C 
C 

WRITE(6,9507) A I F , RE YF , GF 

9507 FORMAT!///, 21X, 'FINAL ESTIMATE... AI =',F12.5,5X, 

* 'REY =' * F8 .1 ,/ ,41X , ' STAB =' , 1PD15.4, ////////) 

GO TO 1 

C 

IOC WRITE( 6,9000 ) 

STOP 

END 
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FUNCTION EIGLS 



PURPOSE 

EIGLS RETURNS THE STABILITY OF THE LEAST STABLE 
EIGENVALUE FUR PLANE POISEULLE FLOW AS DETERMINED 
FRCM A REFERENCE FRAME MOVING IN THE X DIRECTION 
WITH VELOCITY VEL. 

USAGE 

STAB = EIGLS(REY,AR, AI ,BR,BI»VEL) 

DESCRIPTION OF PARAMETERS 

THE FOLLOWING MUST BE SET BY THE CALLING PROGRAM... 
REY,AR,AI ,BR ,BI ,VEL 

REY - THE REYNOLDS NUMBER (REAL*3) 

AR » A I , BR » B I - THE REAL AND IMAGINARY PARTS OF THE 
PERTURBATION WAVE NUMBERS, ALPHA AND BETA 
(REAL*8) 

VEL - THE VELOCITY OF THE COORDINATE SYSTEM FOR 
WHICH THE STABILITY IS DESIRED. (THE MINIMUM, 
AVERAGE, AND MAXIMUM VELOCITIES OF THE FLUID ARE 
O.O, 1.0, AND 1.5 RESPECTIVELY) ( RE AL *8 ) 

THE VALUE OF THE STABILITY OF THE LEAST STABLE 
EIGENVALUE IS RETURNED IN ’EIGLS’. EIGLS MUST BE 
REAL-8 IN THE CALLING PROGRAM. 

OTHER ROUTINES NEEDED 

DEIGEO 



FUNCTION EIGLS(REY,AR,AI ,BR,BI ,VEL) 

IMPLICIT REAL-8 (A-H,0-Z) 

RE AL-3 GR£(60) , G I E ( 60) ,GRO(60) ,GIO (60) 
COMPLEX- 16 XMAT (60 ,30,3) 

CCMF LEX-1 6 YM AT (60,30) ,WVEC(60) , B MAT ( 5,60) 
E GUI VALENCE ( Y MAT ( 1,1) ,XMAT ( 1 , 1 , 3) ) , 

* ( BMAT ( 1 , 1 ) ,XMAT( 1 , 1 * 3 ) ) , 

* (WVECl 1) ,XMAT( 1,6,3) ) 



PC I M = 60 
N = 6C 

CALL SUBROUTINE TO DETERMINE EIGENVALUES FOR BOTH EVEN 
AND ODD EIGENFUNCTIONS. 

CALL CEIGEC(DCMPLX(AR, A I ), DCMPLX(BR,BI ) , REY , N , MD I M, 

* GR E , GI E , GRC , G 1 0 , XMAT , YMAT, BMAT ,WVEC) 

DETERMINE STABILITY FOR EACH EIGENVALUE SAVING THE VALUE 
FOP THE LEAST STABLE ONE. 

TEMP = -1C10 
CO 10 1=1, N 

IF (GRE ( I )+AR*VEL .LT. TEMP) GO TO 5 
TEMP = GRE(I )+AR*VEL 
5 IF ( G PO ( I ) + AR- VEL. LT. TEMP) GO TO 10 
TEMP = GRO(I l+AR-VEL 
10 CONTINUE 

RETURN THE VALUE OF THE MINIMUM STABILITY IN ’EIGLS'. 

EIGLS = TEMP 

RETURN 

END 
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FUNCTION CHM1E1 AND CHM2E1 
(CARTESIAN CCORUINATES) 



FURPCSE 

CHM1E1 AND CHM2E1 RETURN THE VALUES (COMPLEX*16) OF 
ThE COEFFICIENTS FOR THE MATRICES IN THE FINITE 
DIFFERENCE FORM OF THE LINEARIZED NAV I ER-STGKES 
EQUATION FOR POISEUILLE FLOW. BOTH FUNCTIONS RESULT 
FRCM THE LINEAR COMBINATION OF EQUATION 1 AND 
EQUATION 3 TO ELIMINATE THE VELOCITY VECTOR 
POTENTIAL COMPONENT G AND ARBITRARILY SETTING THE 
COMPONENT F TO ZERO. SO, THEY ARE THE COEFFICIENTS 
FOR THE VECTOR POTENTIAL COMPONENT H. CHM2E1 
RETURNS THE TERMS WHICH ARE COEFFICIENTS OF THE 
EIGENVALUE, GAMMA, AND CHM1E1 RETURNS THE REMAINING 
TERMS. 

LSAGE 

XI = CHM IE 1 ( K , Y ) 

X2 = C HM 2E 1 ( K , Y ) 

(CHM1E1 AND CHM2EI MUST BE DECLARED C0MPLEX*16 IN 
CALLING PROGRAM) 

DESCRIPTION OF PARAMETERS 

THE FOLLOWING PARAMETERS MUST BE SET BY THE CALLING 
PR CG RAM 

K - INDICATES THE POINT ON THE FINITE DIFFERENCE 
MESH RELATIVE TO ThE CENTRAL POINT IN THE CENTRAL 
DIFFERENCING SCHEME. IF THE DIFFERENCE IS BEING 
FORMED ABOUT THE N-TH POINT THEN K = 1 REFERS TO 
THE POINT N-2 » K = 2 REFERS TO THE POINT N-l , 

K=3 REFERS TO N, K=4 REFERS TO N + l, AND K=5 
REFERS TO N + 2. 

K - INDICATES WHICH POINT ON THE FINITE DIFFERENCE 
MESH IS REFERRED TO THE CENTRAL POINT. IF THE 
DIFFERENCE IS BEING FORMED ABOUT THE N-TH POINT 
THEN K=1 REFERS TO THE POINT N-2, K=2 REFERS TO 
THE POINT N-l, K=3 REFERS TO N, K=4 REFERS TO N+l, 
AND K=5 REFERS TO N+2. 

Y - THE VALUE OF THE POSITION RELATIVE TO THE CENTER 
OF ThE CHANNEL. THE TWO BOUNDARIES ARE AT Y=+l 
AND Y=- I . 

OTHER ROUTINES NEEDED 
NONE 



FUNCTION C HM 1 E 1 ( K , Y ) 

IMPLICIT COM PL EX#1 6 (A-H,0-Z) 

COMMON / CGEFNT / A, B, G , REY, DEL 
REALMS REY, Y, DEL 

THE FOLLOWING FUNCTIONS (Ml) EVALUATE THE COEFFICIENTS 
OF THE DERIVATIVES OF H FOR ALL TERMS EXCEPT THOSE 
CONTAINING GAMMA. 

CH4MKY) = A/REY 

CH2MKY) = -1 .5DO*A**2*( IDO- Y** 2 ) + 2D0* A* ( A**2 + B**2 ) 

* /REY 

CHOMKY) = -A*( ( a** 2 + B**2)*( i.5D0*A*( lDu-Y**2 )-( A**2 

* +8**2) /REY ) +3D0*A) 
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C THE REMAINING FUNCTIONS ( M2) EVALUATE THE COEFFICIENTS 
C OF THE DERIVATIVES OF H WHICH ARE ALSO COEFFICIENTS 
C CF THE EIGENVALUE GAMMA. 

C 

CF 2M 2 ( Y ) = A 

CH0M2IY) = A*( A**2+B**2) 

SET GF THE FINITE DIFFERENCE VALUES FOR INDEX K FOR MI 

GC TO (1,2, 3, 2,1 ) , K 

1 CHM1E1 = CH4M1(Y)/DEL**4 
GC TO 1GO 

2 CHM1E1 = -4D0*CH4M1(Y)/DEL**4+CH2M1(Y)/DEL**2 
GC TO 100 

3 CHM1E1 = 6C0*CH4M1 (Y )/ DEL**4-2D0*CH2M1 ( Y )/DEL**2 

* + CFOMKY) 

ICC RETURN 

SET LP THE FINITE DIFFERENCE VALUES FOR INDEX K FOR M2 

ENTRY CHM2 El ( K ,Y ) 

GC TO ( 11, 12, 13, 12,11) , K 

11 CHM2E1 = ( CDO, ODO) 

GC TO 200 

12 CHM2E1 = CH2M2(Y)/DEL**2 
GG TO 200 

13 CHM2E1 = -2D0*CH2M2( Y) /DEL**2+CH0M2( Y) 

20C FETUFN 

END 
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SUBROUTINE MULDBM 



FUPFCSE 

MULDBM PERFORMS THE MATRIX MULTIPLICATION BETWEEN A 
SQUAPE MATRIX X AMD A BANDED MATRIX XB WHICH HAS 
BEEN SET UP BY SUBROUTINE BMSET . THE RESULT IS 
PLACED BACK IN X. THAT IS... 

X = (X) (XB) 

LSAGE 

CALL M.ULDBM( X , XB » N » NB, MD I M , TEMPV ) 

CESCRI FT ICN OF PARAMETERS 

THE FCLLOWING MUST BE SET BY THE CALLING PROGRAM. 

N, MDIM, NB,X,XB, TEMPV 

N - TFE SIZE OF THE MATRICES. 

MDIM - THE DIMENSIONS OF THE MATRICES IN THE CALLING 
PRCGRAM. 

N3 - THE NUMBER OF BANDS IN THE BANDED MATRIX, XB. 

X - TFE SQUARE N BY N MATRIX. DIMENSIONED 

( MCI M , MD IM ) IN THE CALLING PROGRAM (COMPLEX* 16 ) 

X9 - TFE N BY N MATRIX WHICH IS STORED IN BANDED 
FORM. DIMENSIONED ( NB , MDIM) IN CALLING PROGRAM. 
(COMPLEX *16) 

TEMPV - A VECTOR WORKSPACE WHICH MUST BE PROVIDED BY 
THE CALLING PROGRAM. DIMENSIONED AT LEAST (N). 
(CCMPL EX*16 ) 

THE FOLLOWING IS OUTPUT BY MULDBM. 

X - SET TO THE PRODUCT OF X AND XB (MDIM, MDIM) 
(CCMPLEX*16) 

RE CU I RED ROUTINES 

NONE 



SUBROUTINE MULDBM ( X, XB ,N ,NB, MDIM ,T EMPV ) 

COMP L EX* lc X ( MD I M, MD IM ) , XB ( NB, MD I M ) , TE MP V< MD I M ) , TEMP 
NBHM = (NB-l)/2 
NBHP = (NB+l)/2 

LCCP CVER INDEX I 

CO 100 1 = 1, N 

STORE COLUMN I OF MATRIX X TEMPORARILY 

DC 10 J=1,N 
10 TEMPV(J) = X ( I , J ) 

FIND PRODUCTS FOR FIRST NBHM SPECIAL CASES, THAT IS WHERE 
WHERE THE BANDED MATRIX DOES NOT HAVE ITS FULL WIDTH 

CC 22 J = 1 , NBHM 
TEMP = ( CCO , ODO) 

JJ = NBFM + J 
DO 21 K= 1 » J J 
JJJ = JJ-K+1 
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21 TEMP = TEMP+TEMPV(K)*X8(JJJ,K) 

22 X ( I , J } = TEMP 

CCMPUTE PRCDUCTS FOR "REGULAR" COMBINATIONS OF ROWS AND 
COLUMNS, THAT IS, THOSE THAT ARE NOT TRUNCATED 
AT TFE BEGINNING OR END BY THE BOUNDARIES 

JF = N-NBHM 
CC 32 J= NBFP , J F 
TEMP = (GCC,0D0) 

CO 31 K = 1 , NB 
JJJ = NB-K+1 

31 TEMP = TEMP4-TEMPVI J-NBHP+K } *XB ( JJJ , J-NBHP+K) 

32 X( I , J ) = TEMP 

FIND PRODUCTS FOR LAST NBHM SPECIAL CASES. 

DC 42 J=1 , NBHM 
TEMP = ( ODC, ODO J 
JJ = NB-J 
CC 41 K = 1 , J J 

41 TEMP = TEMF+TEMPV(N-JJ+K)*XB(NB-K+lt N-JJ+K) 

42 X ( I , N-N6HM+ J } = TEMP 
C 

IOC CONTINUE 
C 

RETURN 

END 
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SUBROUTINE MSET 



PURPOSE 

MSET SETS UP A MATRIX FOR THE FINITE DIFFERENCE 

PROBLEM CF POISEULLE FLOW WITH THE PROPER BOUNDARY 

CONDITIONS FOR THE VELOCITY VECTOR POTENTIAL FOR 

VISCCLS FLOW. 

LS AGE 

CALL MSET(X,N,MDIM» MODE » CFMAT ) 

DESCRIPTION OF PARAMETERS 

THE FOLLOWING MUST BE SET BY THE CALLING PROGRAM 
N,MDIM, MODE, CFMAT 

N - THE SIZE OF THE MATRIX. IF MODE=0, N IS EQUAL 
TO THE NUNBER OF POINTS OF THE FINITE DIFFERENCE 
MESH ACROSS THE CHANNEL NOT INCLUDING THE TWO 
BOUNDARIES. IF MODE=l OR MCDE.= 2, N IS EQUAL TO 
CNE-HALF OF THE NUMBER OF INNER POINTS ACROSS THE 
FULL CHANNEL. IN THIS CASE, THE CHANNEL MUST BE 
DIVIDED INTO AN EVEN NUMBER OF POINTS SO THAT N 
WILL BE AN INTEGER. 

MDIM - ThE COLUMN DIMENSION OF THE MATRIX X. MDIM 
MOST BE .GE. N. 

MODE - IF M0DE=0, THE MATRIX IS SET UP FOR THE FULL 
CHANNEL. IF MODE=l OR MUDE = 2, THE MATRIX IS SET 
UP FCR THE HALF CHANNEL AND THE BOUNDARY 
CONDITIONS ARE SET SUCH THAT THE ODD OR EVEN 
EIGENFUNCTIONS, RESPECTIVELY, ARE SOLVED FOR. 

CFMAT - THE NAME OF A FUNCTION SUBPROGRAM WITH TWO 
PARAMETERS, K AND Y, INDICATING WHICH TERM OF THE 
FINITE DIFFERENCING IS DESIRED, AND THE POSITION 
RELATIVE TO THE CENTER OF THE CHANNEL. MUST BE 
DECLARED EXTERNAL IN THE CALLING PROGRAM. 
(CCMFLEX*16) 

THE FOLLOWING IS OUTPUT BY MSET 

X - THE N BY N MATRIX INTO WHICH THE COEFFICIENTS 
CF ThE FINITE DI FF ER ENC IMG ARE PUT. MUST BE 
DIMENSIONED (MDIM, MDIM) IN THE CALLING PROGRAM. 

( COM PL EX* I 6 ) 

NOTES. . . 

THE EIGENVALUES AND VECTORS OBTAINED BY USING MSET 
TWICE, WITH MODE = 1 AND MODE = 2, ARE THE 
SAME AS USING MSET ONCE WITH M0DE=0, BUT WITH 
N TWICE AS LARGE. 

OTHER ROCTINES NEEDED 

NONE - EXCEPT THE FUNCTION SUBPROGRAM NAME PASSED IN 

THE PARAMETER 'CFMAT*. 



SUBROUTINE MSET (X, N,MD IM , MOD E , C FMA T ) 
REAL*8 REY , Y , DEL Y , DF LOAT 
CCMP L EX* 16 A , B , G 

CCMMCN / CCEFNT / A , B , G , REY , DE LY 

COMPLEX* 16 CFMAT 

COMP LEX* 16 X ( MDI M , MD 1 M ) 
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C COMPUTE GRIC SIZE FOR FINITE DIFFERENCE MESH ACROSS 
C HALF CHANNEL CR FULL CHANNEL 
C 

DELV = 2C0/DFL0AT! N+l ) 

IFIMODE.EC.l .OR. MODE. EG. 2) DELY = 2D0/DFL0AT ( 2*N+ 1 ) 

CHECK IF MATRIX DIMENSIONED LARGE ENOUGH 

IF (N.LE.MDIM) GO TO 1 
WRITE(6,9C00 ) 

900 C FORMAT PO* * * ERROR - ARRAYS NOT DIMENSIONED LARGE', 
* » ENOUGH * * *• ) 



STOP 

ZERO ENTIRE MATRIX 

1 CO 10 1=1, N 
DC 10 J=1,N 
10 X(I , J) = ( CDO , 0D0 ) 

DO SPECIAL CASE AT DISTANCE DELY FROM CHANNEL WALL 
INCLUDING BOUNDARY CCNDITICNS 



Y = 1D0-DELY 

X(l,l) = CFMAT! 3,Y J+CFMAT! 1 , Y) 
XI 1, 2) = CFM AT ( 4 , Y ) 

X ( 1 , 3 ) = CFM AT ( 5 ,Y ) 



DO SPECIAL CASE AT DISTANCE 2*DELY FROM CHANNEL WALL 
INCLUDING ECUNCARY CCNDITICNS 

Y = 1DC-2C0*DELY 
X ( 2 , 1 ) = CFMAT ( 2 , Y ) 

X ( 2 , 2 I = CFMAT (3 , Y ) 

X ( 2 , 3 ) = CFM AT ( 4, Y) 

X ( 2 » 4 ) = CFMAT(5,Y) 

DO ALL REGULAR POINTS IN BETWEEN, THAT IS, THOSE VALUES 
OF Y FCR WHICH ALL 5 FINITE DIFFERENCE GRID POINTS ARE 
INTERIOR TC THE CHANNEL 

IL = N-2 
CO 20 1=3, IL 
K = 1-3 

Y = 1CC-CELY*DFLCAT( I ) 

DO 20 J= 1 , 5 

20 X( I , K + J ) = CFMAT ( J »Y ) 

C 

C FINALLY DO THE TWO SPECIAL CASES WHICH OCCUR EITHER AT 
C THE CENTER OF THE CHANNEL OR AT THE OTHER WALL, DEPENDING 
C CN THE VALUE OF MODE. BOUNDARY CONDITIONS ARE SET UP 
C DEPENDING CN MODE 
C 

Y = 1CC-CELY*DFL0AT(N-1) 

X ( N- 1 , N-3 ) = CFMAT (1 ,Y) 

X ( N- 1 , N-2) = CFMAT { 2 , Y ) 

X ( N- 1 » N- 1 } = CFMAT ( 3 , Y ) 

X { N- 1 , N ) = CFMAT ( 4 , Y J 



IF(MCDE.EQ.l) X ( N- 1 , N ) 
IF ( MCDE. EQ.2 ) X(N-1,N) 

Y = ID C- DELY *D FLOAT! N ) 
X ( N , N- 2 ) = CFMAT! 1 ,Y) 
X! N, N- I ) = CFMAT ( 2 ,Y) 
IF(MCDE.EQ.l) X ( N , N - 1 ) 
IF ! MCDE. EQ.2 ) X(N,N-1) 



CFMAT(4,Y)-CFMAT! 5,Y) 
CFMAT (4,Y)+CFMAT!5,Y) 



CFMAT! 2,Y )-CFMAT( 5,Y) 
CFMAT! 2, Y J + CFMAT! 5,Y) 



X ! N , N J = CFMAT ( 3 , Y ) +CFMAT! 5 , Y) 

IF (MODE. EQ.l ) X ( N » N ) = C FM AT ( 3 , Y ) - CFMAT ( 4 , Y ) 
I F ( MOD E . EG . 2 J X!N,N) = CFMAT l 3 , Y ) + CFMAT ( 4, Y J 

RETURN 

END 
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SUBROUTINE BMSET 



PURPOSE 

THE PURPOSE OF BMSET IS EXACTLY THAT OF MSET EXCEPT 
THAT BMSET TAKES ADVANTAGE OF THE BANDED NATURE OF 
THE FINITE DIFFERENCE MATRICES TO CONSERVE SPACE. 

USAGE 

CALL BMSET ( X , N , MDI M, MODE » CFMAT ) 

DESCRIPTION OF PARAMETERS 

THE PARAMETERS FOR BMSET ARE THE SAME AS THOSE FOR 
MSET WITH THE EXCEPTION THAT THE MATRIX X MUST BE 
DIMENSIONED (5,MDIM) IN THE CALLING PROGRAM. 

NOTE: THE PROCEDURE IS IDENTICAL TO THAT OF MSET. 

COMMENTS HAVE THEREFORE NOT BEEN INCLUDED IN BMSET. 



SUBROUTINE BMSET (X,N,MDIM, MODE, CFMAT ) 

PEAL*8 REY,Y,DELY,DFLOAT 
CCMP L EX* 16 CFMAT 
CCMPLEX-16 A,B,G 
COMPLEX^ 16 XI5.MDIM) 

COMMON / CCEFNT / A , B , G , R EY , DE LY 
CELY = 2DO/D FLOAT! N+ 1 ) 

IF! MODE. EG. I. OR. MODE. EG. 2) DELY = 2DO/DF LOAT { 2*N + 1 ) 

IF (N.LE.MDIM ) GO TO 1 
WRITE ! 6,9000 ) 

9000 FORM AT ! ' 0^ * * ERROR - ARRAYS NOT DIMENSIONED LARGE', 
* • ENOUGH * * *') 

STOP 

1 CO 10 1=1,5 
CO 10 J= 1 , MD IM 
10 X! I , J) = ( CDO , ODO ) 

Y = 1D0-DELY 

X ( 3 , 1 J = CFMAT (3»Y ) + CFMAT ( 1 , Y ) 

X ( 4 , 1 J = CFMAT ( 4, Y) 

X ( 5 , 1 ) = CFMAT t 5 , Y ) 

Y = 1CC-2D0*DELY 

X ! 2 , 2) = CFMAT ( 2 , Y ) 

X ( 3 , 2 ) = CFMAT (3, Y) 

X ( 4 , 2 J = C FM AT { 4 , Y ) 

X ( 5 , 2 ) = CFMAT ( 5 , Y ) 

IL = N-2 
CC 20 1=3, IL 

Y = 1CC-DELY*DFL0AT! I) 

CC 20 J=l,5 

20 X t J , I ) = CFMAT! J,Y) 

Y = 1DC-DELY*DFL0AT( N-l) 

X ( 1 » N- 1 ) = CFMAT { 1 ,Y ) 

X ! 2 , N- 1 ) = CFMAT ! 2 , Y ) 

X ( 3 , N- 1 ) = CFMAT ! 3 , Y J 

X ! 4 , N - 1 J = CFMAT !4 , Y ) 

IF (MODE. EQ. 1 ) X ( 4 , N- 1 ) = CFMAT! 4 ,Y )-CFMAT( 5, Y) 

IF! MG CE. EG. 2) X(4,N-1) = C FMAT ( 4 , Y ) + CFMA T( 5 , Y ) 

Y = 1CC-DELY*DFLGAT( N) 

X! 1 , N ) = CFMAT! 1 ,Y) 

X t 2 , N J = C FM AT ( 2 , Y ) 

IFtMCDE.EC.l ) X t 2 , N ) = C FM A T ( 2 , Y ) -CFMAT ( 5 , Y ) 

I F ! MODE . EQ .2 ) X(2,N) = CFMAT < 2 , Y ) +CFMAT < 5, Y) 

X ! 3 , N ) = CFMAT { 3 ,Y ) +CFMAT ( 5 , Y ) 

IF (MODE. EQ.l ) X ( 3 , N ) = C FM A T ! 3 , Y ) -CF MAT ( 4, Y ) 

IF! MODE. EG. 2) X!3,N) = C FMAT ( 3 , Y J + CFMAl ( 4, Y ) 

RETURN 

END 
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SUBROUTINE DSPLIT 



PURPOSE 

DSPLIT TAKES A MATRIX OF CCMPLEX*16 NUMBERS AND 

SPLITS IT INTO TWO MATRICES, ONE CONTAINING THE REAL 

FART CF THE ORIGINAL MATRIX, AND ONE CONTAINING THE 

IMAGINARY PART. 

USAGE 

CALL CSPLIT (N,MDIM,A,AREAL,AIMAG) 

DESCRIPTION OF PARAMETERS 

N - THE SIZE OF THE MATRIX A, AN N BY N SQUARE 
MATRIX. 

MDIM - THE COLUMN DIMENSION OF MATRIX A 

A - THE INPUT MATRIX. MUST BE DIMENSIONED MDIM BY 
AT LEAST N IN THE CALLING PROGRAM ( COMPLEX* 16 ) 

AREAL, AIMAG - THE OUTPUT MATRICES CONTAINING THE 
REAL AND IMAGINARY PARTS, RESPECTIVELY, OF 
MATRIX A. MUST BE DIMENSIONED (MDIM, MDIM) IN THE 
CALLING PROGRAM. 

NOTES. .. 

MATRIX A AND MATRIX AREAL MAY OVERLAP IF THEY ARE 
DIMENSIONED IN THE CALLING PROGRAM AS FOLLOWS... 

CCMPLEX*16 A( MDIM, MDIM) 

REALMS AREAL ( MDIM, MDIM) , A I M AG ( M DI M , MD I M ) 

EQUIVALENCE ( A ( I , 1 ) , ARE AL( 1 , 1 ) ) 

CTFER ROUTINES NEEDED 

NONE 



SUBROUTINE DSP L I T ( N, MD I M , A , AR , A I ) 

REAL*8 A ( 2, MDIM, MDIM) , AR ( MDIM, MDIM) , AI ( MDI M, MDIM) 
C 

DC 1 J = 1 » N 
CG 1 1=1, N 
A R ( I , J ) = A ( 1 , I , J ) 

1 A I ( I , J ) = A( 2, I , J) 

C 

RETURN 

END 
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SUBROUTINE DEIGEO 



FUFFOSE 

DEIGEC SOLVES THE LINEARIZED N AV IER-ST OKES EQUATION 

FOR POISEULLE FLOW. THE INPUTS TO DEIGEO ARE THE 

WAVE N C S . » ALPHA AND BETA, AND THE REYNOLDS NO. 

DEIGEC OUTPUTS THE EIGENVALUES FOR GAMMA. 

USAGE 

CALL DEIGEO! AL PH A, BETA , REYNG , N ,MDI M, WR EVEN »WI EVEN , 
WRCDD, WIODC,CDM, DM,BM, WV) 

CESCFIPTICN OF PARAMETERS 

THE FOLLOWING MUST BE SET BY THE CALLING PROGRAM... 
ALPHA, BETA, REYNO, N,MDI M 

ALPHA - THE PERTURBATION WAVE NUMBER IN THE FLOW 
DIRECTION (X). (CCMPLEX*16) 

BETA - THE PERTURBATION WAVE NUMBER IN THE DIRECTION 
PERPENDICULAR TO THE FLOW BUT PARALLEL TO THE 
PLATES IZ). ( COMPLEX* 16 ) 

REYNO - THE REYNOLDS NUMBER (REAL*8) 

N - THE SIZE OF THE MATRICES WHICH IS EQUAL TO 
(NC-1I/2 WHERE ND IS THE NUMBER OF DIVISIONS 
ACROSS THE CHANNEL. (NOTE... DEIGEO SOLVES THE 
PROBLEM AC FOSS THE HALF CHANNEL TWICE - ONCE FOR 
THE EIGENVALUES CORRESPONDING TO THE EVEN 
EIGENFUNCTIONS AND ONCE FOR THOSE CORRESPONDING 
TO THE ODD EIGENFUNCTIONS. 

MDIM - ThE COLUMN DIMENSION OF THE MATRICES WHICH 
DEIGEO USES. MDIM MUST BE . GE . N 

THE FOLLOWING ARE OUTPUT BY DEIGEO 
WREVEN ,WI EVEN, WRCDD, W I ODD 

WREVEN »WIEVEN - THE REAL AND IMAGINARY PARTS OF THE 
EIGENVALUES CORRESPONDING TO THE EVEN 
EIGENFUNCTIONS. DIMENSIONED TO AT LEAST N. 
<REAL*8) 



WRCDD, WIODD - THE REAL AND IMAGINARY PARTS OF 
EIGENVALUES CORRESPONDING TO THE ODD 
EIGENFUNCTIONS. DIMENSIONED TO AT LEAST N. 
(REAL*8) 



THE 



THE FOLLOWING 
WORKSPACE. 



MATRICES MUST BE INPUT TO DEIGEO AS 



COM ( MD IM , MDI M) (CCMPLEX*16) 

CM (MDIM, MDIM) (REAL-8) 

BM ( 5 , MD I M) ( CGMPL EX*1 6 ) 

WV(MDIM) ( CCMPLEX-16) 



NOTES. 

THE 

FOR 



MATRICES CAN B 


E OVERLAPPED T 


0 


CONS 


ERVE SPACE, 


EXAMPLE, FOR N 


= 60. 


• • 










CCMPL EX*16 


CDM( 60 


,30 


,3) 








COM PL EX* 1 6 


DM (60, 


30) 


, WV ( 6 


0) 


, BM ( 


5,60) 


EQUIVALENCE 


( DM( 1 


. 1) 


, C D M ( 


1, 


1,3) 


) > 




( BM 1 1 


, 1 ) 


» CD M( 


I, 


1, 3 ) 


), 




W V E C ( 


1) t 


CDM ( I 


,6 


,3) ) 
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NOTE THAT IT IS ONLY THE ACTUAL SIZE OF THESE 
WORKSPACES THAT IS IMPORTANT, NOT THEIR TYPE. 

CTHER FCLTINES NEEDED 

THE FCLLCWING ARE CALLED BY DEIGEO 

CHM1E1»CHM2E2,MSET ,CDMT I N , BMSE T , MULD BM ,DSPLIT, 
EH ESSC jELP.HIC 



SUER CUT I NE DEIGECt AL PH A, BETA, REYNO ,N ,MDI M, 

* WR EVEN , WI E VEN , WR OC D , W I ODD , COM, DM , BM, WV ) 

IMPLICIT CCMPLEX*16( A-H, O-Z) 

DIMENSION IVEC(IOO) 

REAL* 8 WREVEN ( 1) ,WIEVEN( 1) ,WRODD(l ) ,WIODD( 1) 

RE AL*8 C CM ( 1 ) , DM ( 1 ) , BM ( 1 ) , WV ( 1 ) 

RE AL*8 PEY.DELY, REYNO 

COMMON / CCEFNT / A , B , G ,RE Y, DE LY 

EXTERNAL ChMIEl ,CHM2E1 

THIS SUBROUTINE SOLVES THE EQUATION YV = GXV WHERE 
X AND Y ARE MATRICES, V IS THE EIGENVECTOR AND G IS THE 
EIGENVALUE. TEE EIGENVALUES ARE DETERMINED AND PASSED 
BACK TO TEE CALLING PROGRAM IN WRODD, WIODD, WREVEN AND 
WIEVEN. 

A = ALPHA 
E = BETA 
REY = REYNC 

SET UP MATRIX X FOR ODD EIGENVECTORS. 

CALL M SE T( CDM»N , MDI M » 1 ,CHM2E 1 ) 

INVERT MATRIX X. 

CALL CCMTINI N, CDM, MDI M, DETERM) 

SET UP MATRIX Y IN BAND STORAGE MODE FOR ODD EIGENVECTORS. 

CALL BESET (BM,N,MDIM,1 ,CHM1E1) 

MULTIPLY MATRIX Y BY THE INVERSE OF MATRIX X TO CONVERT 
TO TEE STANCARC EIGENVALUE PROBLEM WHICH HAS THE FORM 
(Z-G)V = 0 WHERE Z = < Y) ( I NVERSE l X ) ) . 

CALL MUL CBM ( CDM , BM , N , 5 » MD IK , WV ) 

SPLIT MATRIX INTO REAL AND IMAGINARY PARTS AND CALL 
THE SUBROUTINES TO FIND THE EIGENVALUES. 

CALL D SP L IT( N» MD I M ,CDM ,CDM , DM) 

CALL E HE S S C ( CDM * DM , 1 ,N , N »MD I M , IVEC) 

CALL E LRH1C ( CDM , DM , 1 , N , N , MD I M, WRODD, WIODD, IN ERR, IER) 

IF ( INERR.NE.O) WRI TE ( 6 ,9000 ) INERR , I ER 
900 C FORMAT (' OERROR NUMBER 1 , 17,' ON El GE NV AL UE * , I 7 , / / / ) 

REPEAT TEE SOLUTION FOR EIGENVALUES FOR THE EVEN 
EIGENVECTORS 

CALL MSET(CDM»N,MDIM,2»CHM2E1) 

CALL CD MTIN(N, CDM, MDI M, DETERM) 

CALL BMSET ( BM » N ,MD I M ,2 »CHM I E 1 ) 

CALL MULL'BMl CDM , BM , N , 5 , MD I M , WV ) 

CALL G S P L I T ( N » MD I M , C DM , C DM , D M ) 

CALL EFESSCt CDM, DM ,1 , M , N , MD I M , I V EC ) 

CALL ELRHICt CDM, DM , I ,N ,N ,M DIM, WREVEN, WI EVEN, INERR, IER) 

I F ( INERR.NE.O) WRI TE (6 ,9000) INERR, IER 

RETURN 

END 
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SUBROUTINE CDMTIM (CATEGORY F-l) 
PURPOSE 

INVERT A COMPLEX* 16 MATRIX 
USAGE 

CALL CDMTINtN, A, NDIM, DETERM) 
DESCRIPTION OF PARAMETERS 

N - ORDER OF COMPLEX*16 MA 

(INTEGER) MAXIMUM 'N* 

A - COMP LEX* 1 6 INPUT MATRI 

INVERSE OF ' A' IS RE TU 

NDIM - THE SIZE TO WHICH 'A' 
(ROW DIMENSION OF 'A' 
IN THE DIMENSION S TATM 
CALLING PROGRAM) 

DETEPM - COMPLEX* 16 VALUE OF DE 
RETURNED BY CDMTIN. 



TRIX TO BE INVERTED 
IS 100 

X (DESTROYED). THE 
RNED IN ITS PLACE 

IS DIMENSIONED 
ACTUALLY APPEARING 
ENT OF USER'S 



TERM IN ANT OF 'A' 



REMARKS 

MATRIX 'A' MUST BE A C0MPLEX*8 GENERAL MATRIX 
IF MATRIX 'A' IS SINGULAR THAT MESSAGE IS PRINTED 
' N' MOST BE .LE. NDIM 

SUBROUTINES AND FUNCTIONS REQUIRED 

ONLY BUILT-IN FORTRAN FUNCTIONS 

METHOD 

GAUSSIAN ELIMINATION WITH COLUMN PIVOTING IS USED. 
THE DETERMINANT IS ALSO CALCULATED. A DETERMINANT 
OF ZERO INDICATES THAT MATRIX 'A' IS 
SINGULAR. 



SUBROUTINE CDMTIN ( N , A , NDI M , DETE RM ) 

IMPLICIT REAL*6 (A-H.O-Z) 

COMPLEX* 16 A (NDIM, NO IM) , PIVOT (100) , A MAX, T, SWAP, 
* CE TERM ,U 

I NTE GER*4 IPIVOT( 100) , INDEXl 100,2) 

RE AL*8 TEMP, ALPHA! 100) 

INITIALIZATION 

DETEFM = ( IDO , ODO ) 

DC 20 J= 1 , N 
ALPH A( J) = ODO 
DC 10 1 = 1, N 

10 ALPHA( J)=ALPHA(J)+A( J, I )*DCCNJG( A( J, I ) ) 

ALPHA ( J) =D SORT ( ALPHA (J) ) 

20 IFIVOT (J) =C 
DC 600 1=1, N 

SEARCH FOR PIVOT ELEMENT 



AMAX = 


(ODO, ODO) 








CO 


1 C 5 


J = 1 » N 










IF 


( IP 


I V C T ( J ) ~ 


1 ) 


60, 


105, 


60 


CO 


100 


K= 1 , N 










IF 


( IP 


IVOT(K)- 


1) 


80, 


100, 


740 
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8C TEMP=ANAX*DCON JG(AMAX)-A( J , K ) *DCON JG < A { J ,K)) 
IF(TEMP) £5, 85, 100 

85 IRCW = J 
ICCLUM =K 
AM AX = A ( J , K ) 

ICC CONTINUE 
105 CCNTINUE 

IPIVCT ( I CGLUM) = I PIVOT! ICOLUM )+l 

INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

IF ( IRCW-ICCLUM) 140, 260, 140 

140 CETERM =DETEP.M 
DC 200 L = 1 , N 
SW A F = A < I FCW, L) 

A ( IROW » L ) =A ( ICOLUM , L ) 

20C A ( ICC LLN , L )=SWAP 
SNAP=ALPHA ( I ROW) 

ALPHA ( IR 0 W ) = AL PH A ( ICOLUM) 

ALPHA! ICCLUM ) =SWAP 
260 INCEX < 1 , 1 ) = 1 ROW 

INDEX! 1,2) =1 COLUM 
PIVOT! I) =A( ICOLUM, ICOLUM) 

L= P I VOT ( I) 

TEMP=PIVCT ( I ) *DCCN JG ( P I VOT ! I ) ) 

IF(TEMP) 330, 720, 330 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

330 A! ICOLUM, ICOLUM) = (1D0.0D0) 

CC 350 L=1,N 
L=FI VOT! I ) 

350 A! ICOLUM, L)=A( ICOLUM, L)/U 

REDUCE NCN-PIVOT ROWS 

38C CC 550 L 1 = 1 » N 

IF(Ll-ICCLUM) 400, 550, 400 
400 T = A! LI , ICCLUM) 

A ( Li , I CCLUM) = ! ODO, ODO) 

CC 450 L = 1 , N 
U = A l ICCLUM,L) 

450 A(L1,L)=A!L1,L)-U*T 
55C CCNTINLE 
6CC CCNTINUE 

INTERCHANGE COLUMNS 

6 2 C DC 7 1 C 1 = 1, N 
L=N+1-I 

IF! INDEX(L,1 )-INDEX!L, 2) ) 630, 710, 630 
63 C JROW = I NDEX (L , 1 ) 

JCCLUM=INDEX ( L , 2 ) 

DO 705 K = 1 , N 
SWAP = A(K, JPOW) 

A ( K , JPCW) = A( K, J CCLUM ) 

A(K, JCCLUM)=SWAP 
7C5 CONTINUE 
710 CCNTINUE 
RE TU F N 

720 WRITE(6,73C) 

730 FORMAT ( 20H MATRIX IS SINGULAR) 

74C FETUFN 
END 
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SUBROUTINE DFIT2 



PURPOSE 

THIS SUBROUTINE TAKES 3 POINTS (XI, Y), (X2,Y2), 

AND (X3,Y3), AND USES A PARABOLIC FIT TO FIND THE 
EXTREMAL VALUE OF Y AND THE CORRESPONDING VALUE FOR 

X WHICH ARE RETURNED IN XM AND YM. IT ALSO TAKES AN 
INPUT VALUE, YVAL, AND COMPUTES THE TWO SOLUTIONS 
FOR X TO THE PARABOLIC I NTERPC LA TI ON , TAKES THE 
SOLUTICN NEAREST TO XI, AND RETURNS IT IN XNEXT . 

IF THERE IS NO REAL SOLUTIGN TO THE EQUATION, THEN 
XNEXT IS SET TO 0.0. 

USAGE 

CALL CFIT2 <X1,X2 ,X3,Y1 ,Y 2, Y3, YVAL, XNEXT, XM,YM) 

DESCRIPTION OF PARAMETERS 

THE FOLLOWING MUST BE SET BY THE CALLING PROGRAM 
X1,X2,X3,Y1,Y2,Y3,YVAL 

XI ,X2 , X3 , - THE X VALUES OF 3 POINTS (REAL*8) 
Yi,Y2,Y3 - THE Y VALUES OF 3 POINTS (REAL* 8) 

YVAL - A VALUE OF Y FOR WHICH THE CORRESPONDING 

X VALUE IS DESIRED ( REAL*8 ) 

THE FOLLOWING ARE SET BY DFIT2 
XNEXT, XM,YM 

XNEXT - SET TO THE VALUE OF X NEAREST XI 
CORRESPONDING TO YVAL (REALMS) 

XM , YM - SET TO THE MIN Y VALUE FROM THE SECOND 
ORDER INTERPOLATION, WITH THE CORRESPONDING X 
VALUE ( RE AL* 8 ) 



SUBROUTINE DFIT2(X1,X2,X3,Y1,Y2,Y3 , YVAL , XNEXT , XM , YM } 
IMPLICIT R EAL*8 ( A— Z ) 

A = ( (Yl-Y2)/(Xl-X2)-( Y1-Y3)/(XI-X3) ) 

* /( (Xl**2-X2**2)/(Xl-X2)-(Xi**2-X3**2 ) / ( X 1— X3 ) ) 

B = (Y1-Y2-A*(X1**2-X2**2I ) / (X1-X2) 

C = Y1-A*X1**2-8*X1 
XM = -B/(2D0*A) 

YM = A*XM**2+B*XM+C 
D = B**2-4D0*A*< C-YVAL) 

IF(D.GE.CDC) GO TO 10 
XNEXT = ODO 
RETURN 
C 

1C XN1 = (-B+DSQRT ( D) )/(2D0*A) 

XN2 = (-B-DSQRT(D) )/(2D0*A) 

XNEXT = XN1 

IF(DABS(XN2-X1) .LT.DABS(XNl-Xl) ) XNEXT = XN2 

RETURN 

END 
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PROGRAM #R 1 A 

PROGRAM TO PRINT EIGENVALUES 
FOR THE 3-D CYLINDRICAL FLOW PROBLEM 
N = 0 FOR G 



INPUT VALUES ARE N, REY, AR, AND AI. 
(USING NAMELIST 'LIST* ) 



IMPLICIT REAL*8( A-H, 0- Z ) 

COMPLEX 5 * 16 A, B, DETERM, G 

COMPLEX 5 * 16 XM A T (31,31) ,YMAT(31,31) ,WV(31) 

COMPLEX 5 * 16 CGM1E2 ,CGM2E2 
REALMS G R ( 3 1 ) , G I ( 3 1 ) 

REAL-4 GR4 ( 3 1 ) ,GI4(31) 

DIMENSION IVEC(IOO) 

COMMON / CCEFNT / A , B, G , R EY , DELR 
NAMELIST / LIST / N , REY , AR , A I 
EXTERNAL CGM1E2, CGM2 E2 
C 

MCIM = 31 
N = 30 

REY = 6C00C0 
AR = OCO 
AI = IDO 
B = (OCC.ODO) 

C 

1 READ(5,LIST,END=I00) 

A = DC MPLX (AR , AI ) 

W R I T E ( 6 , 9000 ) 

900C FORMAT Cl') 

WRITE! 6,9001) N , RE Y , A 

9C01 FORMAT ( e COMPONENT G IN EQUATION 2',/, 'ON =',I4,/, 

* • REY = ' , F10» 2, 8X , ' ALPHA =• , 2F1 2. 7,8X , • NBI = 0 s ) 

C 

KP1 = N+l 

CALL SET2(XMAT ,NPI ,1 , MDI M, CGM2 E2 , 1 ,1) 

CALL CCMT IN( N , XMAT , MDI M, DETERM ) 

CALL SET2IYMAT ,NP1 , 1 , MDI M, CGM1 E2 , 1 ,1) 

CALL MULM( XMAT , YMAT, N , MDI M,WV) 

CALL D SP L IT( N , MD I M ,XMAT , XMAT , YMAT) 

CALL EF ESSC( XM AT , YMA T , 1 , N , N , MD I M , I VEC ) 

CALL ELRHiC(XMAT ,YMAT, 1 ,N,N,MDI M,GR,GI , INERR, IER) 
IF( INERR. NE.O) WRI TE ( 6 , 901 0 ) INERR, IER 
9010 FORMAT ( ' 0 5 * 5 * * ERROR NUMBER', 17,' ON EIGENVALUE', 

=* 1 7, • * * *' ,///) 

C 

WRITE(6,9CC2) 

900 2 FORM AT (///,16X,' GAMMA REAL ', 5X ,' GAMMA I MAG') 

DC 10 1=1, N 

GR4 ( I ) = SNGL ( GR ( I ) ) 

G 1 4( I) = SNGL ( G I ( I ) ) 

10 WR I T E ( 6 , 9003 ) I , GR ( I ) , G I ( I ) 

9003 FORMAT ( 'O' , 1 10 , 1 F2D1 5. 4) 

C 

WR IT E ( 6 ,9000 ) 

CALL F LOT P ( GR4 ,G 14 ,N , 0 ) 

WRITE(6,9G01 ) N , RE Y , A 
C 

GO TO 1 

IOC WR IT E ( 6 , 9000 ) 

STOP 

END 
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PROGRAM #R1B 

PROGRAM TO PRINT EIGENVALUES 
FOR THE 3-D CYLINDRICAL FLOW PROBLEM 
N = 0 FOR H 



INPUT VALUES ARE N, REY, AR, AND AI. 
(USING NAMELIST 'LIST' ) 



IMPLICIT REAL*8( A-H,C-Z) 

COMPLEX* 16 A, B, DETERM, G 

COMPLEX* 16 XMAT(31 ,31 ) ,YMAT(31»31) ,WV(31) 

COMPLEX* 16 CHM1E1 ,CHM2E1 
RE AL*8 GR<31) ,GI (31) 

RE AL*4 GR4 ( 3 i ) ,GI4(31) 

DIMENSION IVEC(IOO) 

COMMON / CCEFNT / A , B , G , RE Y , DE LR 
NAMELIST / LIST / N , REY , AR , A I 
EXTERNAL CHM1 El , CHM2 E 1 
C 

MDIM = 31 
N = 30 

REY = 6COODO 
AR = OCO 
AI = IDO 
B = (OCO, OCO) 

C 

1 READ (5,LIST,END=100) 

A = DCMPLX (AR, AI ) 

WRITE(6,9COO) 

9000 FORMAT Cl*) 

WRITE (6 , 9001 ) N , REY , A 

9001 FORMAT ( 1 COMPONENT H IN EQUATION l',/,-'ON =',I4,/, 

* ' REY= ' »F10.2,8X,'ALPHA = • , 2F1 2. 7 , 8X , • N3 1 = O') 

NP1 = N+l 

CALL S ET2( XMAT »NP 1 ,1 , MDI M, CHM2E 1 ,1 ,1) 

CALL CCMT INI N, XMAT , MCI M, DETERM) 

CALL SET2( YMAT ,NP1 , 1 , MDI M, CHM1 E 1 , I , 1 ) 

YMAT (N ,N) = YMAT ( N ,N ) + CHMI El ( 5 ,DELR) 

CALL M U L M ( XM A T , Y M A T , N , M D I M , W V ) 

CALL DSP L IT( N, MDIM ,XMAT, XMAT , YMAT) 

CALL EFESSC( XMAT, YMAT, 1 , N , N , MD I M , I VE C ) 

CALL ELRH 1C (XMAT, YMAT, 1 , N, N , MD I M ,G R , G I , I NE RR , I ER ) 

IF ( INERR.NE. 0) W PI TE ( 6 ,901 0) INERR, I ER 
9010 FORM AT ( ' 0* * * ERROR NUMBER', 17,' ON EIGENVALUE', 

* 17, ' * * *',///) 

C 

WRITE (6, 9002) 

9002 FCRMAT(///,16X, 'GAMMA REAL 5X GAMMA I MAG' ) 

CC 10 1=1, N 

GR4 ( I) = SNGL ( GR ( I )) 

GI4( I) = SNGL ( G I (I ) ) 

10 WRITE(6,9003) I , GR ( I ) , GI ( I ) 

9003 FORMAT ( 'O' ,1 10 , 1P2D15.4) 

C 

WRITE( 6,9000) 

CALL FLOTP(GR4,GI4,N,0) 

WRI T E ( 6 , 9001 ) N, REY , A 
C 

GO TO 1 

IOC KR IT E ( 6 , 9000 ) 

STOP 

END 
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PROGRAM #R2 

PROGRAM TO PRINT EIGENVALUES 
FOR THE 3-D CYLINDRICAL FLOW PROBLEM 

N = 1 



IMPLICIT REAL*8( A-H, C-Z) 

CCMFLEX*16 A, B, DETERM, G 

COMPLEX* 16 XM AT (60,60) ,YMAT{ 60,60) ,WV(60) 

COMPLEX* 16 CFM1E1,CFM2E1 ,CFM 1 E2 , CF M2E2 ,CGM1E1 ,CGM2E1 , 

* CGM1E2,CGM2E2,CHM1E1 ,CHM2EI , CriMl E2 ,CHM2E2 
RE AL*8 GR(60) ,GI(6C) 

RE Al*4 G R4 ( 60 ) ,GI4(60) 

Cl ME NS ION IVEC(IOO) 

COMMON / CCEFNT / A , B , G , REY , DELR 
NAMELIST / LIST / N , REY , AR , A I 

EXTERNAL CFM 1E1 , CFM2 El , C FM1 E2 , CFM2 E2 ,C GM I E 1 ,CGM2E I , 

* CGM1E2,CGM2E2,CHM1 El ,CHM2E1 , CHM1 E2 ,C HM2E2 

MDIM = 6C 
N = 30 

FEY = 6C00D0 
AR = OCO 
A I = ICO 
B = ( OCC , IDO ) 

C 

1 READ ( 5,LIST,END=100) 

A = CCMPLX (AR, AI ) 

WRITE ( 6 , SOOO ) 

90CC FORMAT Cl') 

WR IT E ( 6 , 9001 ) N , REY , A 

9001 FORMAT ( ' N =',I4,/,' REY = • , F 10 . 2 , 3X , * AL PH A = ',2F12.7, 

* 3X,'NBI = 1',/,' FOR VELOCITY VECTOR POTENTIAL ', 

* »H SET EQUAL TC ZERO' ) 

C 

MS I ZE = 2*<N-1) 

CALL SET2(XMAT,N,1 , MDI M , CFM2 El , 1,1 ) 

CALL SET2(XMAT ,N,1 ,MDI N,CGM2E1 , 1,N) 

CALL S ET 2 ( XMAT , N , 1 ,MD I M , CFM2 E2 , N , 1 ) 

CALL S ET 2 ( XMAT , N , 1 , MC I M , CGM2 E2 , N ,N ) 

C 

XMAT (N-l ,N-1 ) = XMAT ( N-l , N- I ) + CFM2 El ( 5 , DELR ) 

XM.AT (M SIZE, N-l ) = XMAT ( MS I ZE , N- 1 ) + CFM2 £2 ( 5 , DELR ) 

C 

CALL CDMTINt MSI ZE, XMAT, MDIM, DETERM ) 

C 

CALL SET2 ( YMAT ,N , 1 ,MDI M, CFM1 El , 1 , 1 ) 

CALL SET 2 ( YMAT ,N , 1 , MDI M , CGM 1 E 1 , i ,N) 

CALL SET2 ( YMAT ,N , 1 , M D I M , CFM1 E2 , N , 1 ) 

CALL SET2( YMAT, N,1 , MDIM, CG Ml E2,N,N) 

C 

YMAT (N-l , N-l ) = YMAT ( N-l , N-l ) + CFM1 El ( 5 ,DELR) 

YMAT ( M S I ZE »N— 1 ) = YMAT ( MSI ZE , N- 1 ) + CF Ml E2 ( 5 , DEL R) 

C 

CALL MULM ( XMAT , YMAT, MS I Z E, MD IM , WV ) 

CALL D SPL I T( MS I ZE , MD I M ,X MAT, XMAT , YMAT) 

CALL ERE SSC( XMAT , YMA T , 1 , MS I Z E , M S I Z E , MDI M , I VEC ) 

CALL ELRH1C( XMAT, YMAT, 1, MSI ZE, MSI ZE, MDIM, 

* GR,GI,INERR, IER) 

I F ( INERR.NE.O) WRI TE ( 6 ,901 0 ) I NERR , I ER 
9010 FORMAT ( • 0* * * ERROR NUMBER', 17,' ON EIGENVALUE', 

* 17,' * * * ' ,///) 

L WRITE! 6,9002 ) 

9002 FORMAT (///,16X, 'GAMMA REAL * , 10X , 'G AMMA IMAG') 

DC 10 I = 1 , MS I Z E 

GR4CI) = SNGL ( GR ( I )) 

G I4( I) = SNGL ( G I ( I ) ) 

10 WR IT E ( 6 , 9003 ) I , GR ( I ) , GI ( I ) 
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9003 FORMAT ( ' C * »I 10ilP2D2C. 10) 

C 

VnRITE(6,9000) 

CALL PLGTR(GR4,GI4,MSIZE,0) 
ViRITE ( 6 > 9001 } N,REY,A 
C 

GC TO 1 

IOC HR IT E ( 6 » 9000 ) 

STOP 

ENO 
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FUNCTICN CHM1E1(K,R) 

IMPLICIT CCMPLEX*16 (A-H,0-Z) 

CCMfICN / CCEFNT / A , B , G , REY , DEL 
REALMS REY, R, DEL 
C 

T1(R) = A**2+( B/R) **2 
T 2 (R ) = T1(R) /REY-A*2D0*( 1D0-R**2) 

C 

CF4MHR) = B/( R*REY) 

CF3MKR) = 2D0*B/( R**2*REY) 

CF2MKR) = 8*T2(R)/R-6/(R**3*REY)+B*Tl(R)/(R*REY) 
CF1MKR) = B*T2l R) /R**2+ 6/-V R**4* RE Y ) -3D0*B**3/ ( R**4 

* *REY) +A**2*B/ <R**2*REY) 

CFOMl(R) = T2 ( R ) *8 5t T 1 ( R ) /R + 400*6**3/ ( R** 5*RE Y ) 

C 

CF2M2 ( R ) = B/R 
CF 1M2 ( R ) = B/R** 2 
CF0M2 ( F ) = B*T1 ( R) /R 
C 

CG2MKR) = -4D0*A*B/ (R**2*REY) 

CG1MKR) = 4D0*A*B/( P**3*REY ) 

CGOMl(R) = -4D0*A*3/(R**4*REY)-2D0*A*8*(T2(R) 

* +T1(R)/REY)/R**2 
C 

CG0M2( R ) = -2D0*A*B/R**2 
C 

CH4M1 ( P) = -A/REY 
CH3MKP) = -2 DO* A / ( R*R EY J 

CH2MKR) = - A*T1 ( R ) / REY- A*T2 ( R ) + 3D0*A/ ( R**2*REY ) 
CH1MKR) = 3D0*A*B** 2 / { R**3*REY ) -A*T 2 ( R) /R-3D0*A 

* / ( R**3*REY ) -A* *3/ ( R*REY ) 

CHOMl(R) = T2(R)*( A/R**2-A*T1( R) )+A*B**2/{ R**4*REY) 

* +3L0*A/(R**4*REY )+A**3/(R**2*REY) 

C 

CH2M2 ( R ) = -A 

CF 1M 2 ( R ) = -A/R 

CH0M21P) = A/R**2-A*T1 (R ) 

C 

GO TG (11,12, 13, 14,15) ,K 

11 CHM1E1 = CH4M1(R)/DEL**4-CH3M1 ( R ) / ( 2 DO *D EL** 3) 

GC TC ICO 

12 CHM1E1 = -4D0*CH4M 1 ( R. ) /DEL**4+ 2D0*CH3M 1( R) / ( 2 DO* DEL** 3 

* ) +CH2M1 (R)/DEL**2-CH1M1 (R)/(2D0*DEL) 

GC TC ICO 

13 CHM1E1 = 6D0*CH4M1(R)/DEL**4-2D0*CH2M1(R)/DEL**2 

* + CFOPKR) 

GO TC IOC 

14 CHM1E1 = -4D0*CH4M 1( R ) /DE L**4- 2D0* CH3M1 ( R) / ( 2D0*DEL**3 

* ) +CH2K1 (R)/DEL**2+CH1M1 (R)/( 2D0*DEL) 

GO TC ICO 

15 CHM1E1 = CH4M 1 ( R ) /DE L**4+CH3M 1 ( R ) / ( 2 DO*D EL**3 ) 

IOC RETURN 

C 

ENTRY CHM2E1 ( K,R ) 

GC TC (21,22,23,24,21) ,K 

21 CFiMZEl = (ODO, ODO) 

GC TO 200 

22 CHM2E1 = CH2M2(R)/DEL**2-CH1M2 (R)/ (2D0*DEL) 

GO TC 200 

22 CHM2E1 = -2D0*CH2M2(R) /DEL**2+CH0M2( R) 

GC TO 200 

24 CHM2E1 = CH2M2(R)/DEL**2+CH1M2(R)/(2D0*DEL) 

200 RETURN 
C 

ENTRY CGM1E1 ( K , R ) 

GC TC (31,22,33,34,31) ,K 

31 CGK1E1 = (ODO,ODO) 

GC TC 300 

32 CGK1E1 = CG2M1 (R)/DEL**2-CG1M1 ( R )/ ( 2D0*DEL ) 

GC TC 3C0 

33 CGM1E1 = - 2D0*CG 2M 1 1 R) /DEL** 2+CG OM 1 ( R) 

GC TC 300 
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34 CGM1E1 = CG2M1(R)/DEL**2+CG1M1(R)/ (2D0*DEL) 

3CC RETURN 

ENTRY CGM2E1(K,R) 

GC TG (41 ,41 ,42 ,41 ,41 ) ,K 

41 CGM2 E 1 = (000,000) 

GO TO 400 

42 CGM2 El = CG0M2 ( R ) 

4CC RETURN 

ENTRY CFR1E1(K,R) 

GC TO (5 1, 52,53, 54 ,55) ,K 

51 CFM1E1 = CF4M1(R)/DEL**4-CF3M1(R)/(2D0*DEL**31 
GC TC 500 

52 CFM1E1 = -4D0-CF4M1 ( R) /D EL**4+2D0*CF3M 1 ( R) /( 2DO*DEL**3 

* ) +CF2R1(R)/DEL**2-CF1M1 (R)/(2D0*DEL) 

GC TC 500 

53 CFM1E1 = 6C0*CF4M1 ( R ) /DEL**4-2D0*C F2M1 (R )/DEL**2 

* +CF0MKR) 

C-C TC 500 

54 CFM1E1 = -4D0*CF4M1(R)/DEL**4-2D0*CF3MMR) /(2D0*DEL**3 

* ) +CF2M1 (R)/DEL**2+CF1M1 ( R )/ ( 2D0*DEL ) 

G C T C c C 0 

55 CFMlEl"= CF4M1(R)/DEL**4+CF3M1(R)/(2D0*DEL**3) 

500 RETURN 

ENTRY CFM2E1(K,R) 

C-C TC (61, 62, 63, 64, 61 ), K 

61 CFM2E1 = (000,000) 

GC TO 600 

62 CFM2E1 = CF2M2(R)/DEL**2-CF1M2(R)/(2D0*DEL) 

GC TO 600 

63 CFM2E1 = -2D0*CF2M2(R)/DEL**2+CF0M2(R) 

GC TC 600 

64 CFM2E1 = CF2N2(R)/DEL**2+CF1M2(R)/ (2D0*DEL) 

60 C RETURN 

END 
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c 

c 
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FUNCTION C HM 1 E 2 ( K » R ) 

IMPLICIT CCMPL£X*16(A-H,0-Z) 

CCMMCN / CCEFNT / A , B , G , REY , DE L 
REAL*8 REY ,R ,DEL 

T1(R) = A**2+( 3/R)**2 

T2(R) = T1(R)/REY-A*2D0*(1D0-R**2) 

CF3MKR) = A/REY 
CF2MKR) = A/(R*REY) 

CF1MKR) = A*( T2 (R)-1D0/ (R**2*REY) ) 

CFOMl(R) = 4DO*B**2/R-2DO*A*6**2/{ R**3*REY ) 

CF1M2 ( R ) = A 

CG2M1( R) = -TKRJ/REY 

CG1MKF) = 2D0*B**2/(R**3*REY)-T1( R)/(R*REY) 
CGOKl(R) = T1(R)*(1D0/(R**2*REY)-T2(R) )-2DO*B**2 

* / ( R*$4*REY ) 

CG0M2( R) = — T 1 ( R ) 

CH3MKR) = B/ ( R*REY ) 

CH2MHR) = 2D0*B/( R**2*REY) 

CH1MKR) = B*(T2(R)-1D0/(R**2*REY) )/R 
CHOMl(R) = -4D0*A*B+B*T2(R)/R**2+B/(R**4*REY) 

* +2CC*A**2*B/ ( R**2*REY) 

CH1M2 ( P ) = B/R 
CHOM2(R) = B/R** 2 

GC TC (11, 12, 13, 14, 15), K 

11 CHM1E2 = -CH3MK R) /( 2D0*DEL**3) 

GC TO 100 

12 CHM1E2 = 2D0*CH3Ml(R)/( 2 DO “DEL** 3) +CH2M1 (R)/DEL**2 

* -CH1M1(R)/(2DO*OEL) 

GC TO 100 

13 CHM1E2 = -2D0*CH2M1(R)/DEL**2+CH0M1(R) 

GO TO 100 

14 CHK1E2 = -200*CH3Ml(R)/( 2D0*CEL**3 )*CH2M1( R) /DEL**2 

* *CH1M1(R)/(2D0*DEL) 

GC TO ICC 

15 CHM1E2 = CH3M1(R)/(2C0*DEL**3) 

IOC RETURN 

ENTRY CHM2 E2 ( K , R ) 

GC TO (21,22,23,24,21) ,K 

21 CHM2E2 = (000,000) 

GC TC 2C0 

22 CHM2E2 = -CH 1M2 ( R) / ( 2D0*DE L) 

GC TO 200 

23 CF.M2E2 = CH0M21R) 

GC TC 200 

24 CHM2E2 = CH1M2 ( R ) / ( 2D0*DEL ) 

200 RETURN 

ENTRY CGM 1 E2 ( K ,P ) 

GC TO (31,32,33,34,31) ,K 

31 CGM1E2 = (000,000) 

GC TC 300 • 

32 CGM1E2 = CG2M1 (R)/DEL**2-CG1M1 ( R ) / ( 200*0 EL ) 

GO TO 3C0 

33 CGM1E2 = -2D0*CG2M1(R)/DEL**2+CG0M1(R) 

GC TO 3CC 

34 CGM 1 E 2 = CG2M1(R)/DEL**2+CG1M1(R)/(2D0*DEL) 

30C RETURN 

ENTRY CGM2E2 ( K , R) 

GC TO (4 1,41 ,42,41 ,41 ) ,K 

41 CGM2E2 = (000,000) 

GC TO 4CC 

42 CGM2E2 = CG0M2(R) 
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40C RETURN 



ENTRY CFM1E2 ( K ,R ) 

GC TO (5 1, 52,53, 54,55) ,K 

51 CFRIE2 = -CF3M1 (R)/{ 2DO*DEL**3 ) 

GG TG 500 

52 CFM1E2 = 2D0*CF3M1(R)/(2D0*DEL**3) +CF2M1 (R)/DEL**2 

* -CF1M1(R)/(2D0*DEL ) 

GO TO 5C0 

53 CFM1E2 = -2D0*CF2M1(R) /D EL**2 + CF OM 1 ( R) 

GC TG 500 

54 CFM1E2 = -2D0*CF3M1(R)/(2D0*DEL**3)+CF2M1(R)/DEL**2 

* +CFlf'l(R)/(2C0*0EL) 

GC TO 5CC 

55 CFM1E2 = CF3M1(R)/(2D0*DEL**3> 

50C RETURN 

ENTRY CFM2 E2 ( K » R ) 

GC TC (61,62,61, 64 ,61) ,K 

61 CFM2E2 = ( ODO , ODO ) 

GG TO 600 

62 CFR2E2 = -CF 1M2 ( R ) / ( 2D0* DE L ) 

GC TG 600 

64 CFM2E2 = CF1M2(R)/(2D0*DEL) 

600 RETURN 
END 
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SUBf OUTINE SET2I X,N, MI NR , MDI M , CF MA T , NQV , NQH ) 

RE AL*8 REY,R,DEL,DFLOAT 

CCMPLEX* 16 A , 8 , G 

COMMON / CCEFNT / A , B , G , RE Y , DE L 

CCMPLEX*16 CFMAT 

CC M PLE X* 1 6 X(MDIMjMUIM) 

CEL = 1CO/DFLOAT ( N) 

MV = NCV-1 
MH = NQH- I 

IF(N+MV.LE.MDIM.AND.N+MH.LE.MDIM) GO TO 1 
to RI T E ( 6 , 9000 ) 

9000 FORMAT ( ' 0* * * ERROR - ARRAYS NOT DIMENSIONED LARGE', 
* ' ENOUGH * * *') 

STOP 

1 CC 10 1=1, N 
CO 10 J= 1 , N 

10 X( I + MV ,J +MH ) = ( ODO, ODO) 

R = 1C0-CEL 

X ( 1 + MV , 1 +MHJ = CFMAT (3 , R } +CFMAT ( 1, R) 

X ( 1+MV , 2 +MH) = CFMAT (4 ,R ) 

X ( 1+ M V , 3 + MH) = C FM AT ( 5 , R ) 

R = 1 D 0- 2 D 0* D E L 
X ( 2+ M V , 1 +MH ) = CFMAT ( 2 , R ) 

X ( 2+ MV , 2 + MH) = C F M A T ( 3 , R ) 

X { 2+MV , 3 + MH) = C FM AT ( 4 , R ) 

X ( 2+MV , 4+Mh) = C FM AT ( 5 , R ) 

IL = N-2 
CO 20 1 = 3, IL 
K = 1-3 

R = 1D0-DEL-DFL0ATI I ) 

DO 20 J = 1 , 5 

20 X( I+MV, K+J+MH) = CFMAT ( J , R ) 

R = DEL 

X ( N- 1+MV ,N-3+MH ) = CFMAT(l.R) 

X(N-l + MV,N-2+MH) = C FM AT ( 2 , R ) 

X ( N- 1 + MV , N— 1 + MH ) = CFMAT ( 3 , R ) 

X( N- 1+MV , N+MH) = CFMAT ( 4 , R ) 

IF(MINP.NE.O) GO TO 30 
F = CDC 

X ( N+M V »N-2+MH ) = CFMAT (1,R) 

X(N+MV,N-1+MH) = C FM AT ( 2 , R ) 

Xt N+MV.N+MH) = CFMAT ( 3 , R) 

30 RETURN 
END 
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